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STATUS  OF  THE  RESEARCH  EFFORT 

The  research  reported  here  is  in  the  area  of  numerical  analysis. 
There  are  four  parts  to  this  report: 

1.  Computational  error  estimates  and  deferred  corrections 
for  differential  and  Integral  equations 

II.  Roundoff  error  for  variants  of  Gaussian  elimination 

III.  Multistep  methods  for  ordinary  differential  equations 

IV.  Multi-grid  methods  for  elliptic  partial  differential 
equations 

A  more  detailed  outline  is  given  by  the  table  of  contents,  which  follows. 

The  work  described  in  this  report  falls  into  three  categories: 
(i)  work  that  is  reported  elsewhere,  for  which  we  include  little  more  than 
an  abstract  and  an  outline,  (ii)  work  that  will  not  be  reported  elsewhere, 
for  which  much  more  detail  is  provided,  (iii)  work  that  is  not  yet  ready 
for  publication,  for  which  unpolished,  incomplete  results  are  given. 
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I.  COMPUTATIONAL  ERROR  ESTIMATES  AND  DEFERRED  CORRECTIONS  FOR 
OPERATOR  EQUATIONS 

Digital  computer  simulations  of  physical  phenomena  often  involve 
the  numerical  solution  of  partial  differential  equations  with  initial  and 
boundary  conditions.  It  is  clearly  Important  to  have  an  estimate  of  the 
accuracy  of  the  computed  solutions.  We  have  identified  ten  ways  to 
estimate  the  error  in  the  numerical  solution,  often  called  the  global 
error.  Good  theoretical  and  practical  results  for  estimating  error  for 
smooth  problems  were  obtained  by  Lindberg  (1976)  using  a  generalization 
of  Fox's  deferred  difference  correction.  The  idea  is  to  produce  another 
numerical  solution  of  enhanced  accuracy  by  subtracting  out  from  the 
difference  equation  of  the  improved  solution  local  error  estimates 
computed  for  the  original  solution.  We  have  since  been  able  to  simplify 
and  strengthen  Lindberg 's  theorems.  The  result  is  a  theoretical  framework 
for  proving  accuracy  results  for  deferred  correction.  Application  of 
this  theory  to  ordinary  differential  equations  yields  useful  theoretical 
results.  The  practical  benefit  is  to  offer  guidance  in  the  construction 
of  local  error  estimators  by  identifying  those  properties  that  are 
Important  for  accurate  error  estimation.  These  ideas  have  been  applied  to 
time-dependent  partial  differential  equations  with  nonsmooth  solutions, 
particularly  hyperbolic  problems  with  shock  discontinuities.  For  these 
problems  it  would  be  very  difficult  to  obtain  error  estimates  that  can 
be  theoretically  justified;  we  were  content  to  devise  and  test  schemes 
having  only  limited  theoretical  support. 

An  introduction  to  the  idea  of  deferred  correction  and  our 
style  of  error  analysis  is  given  in  Skeel  (1980,  section  1). 
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1.  Error  Estimation  and  Iterative  Improvement  for  the  Numerical  Solution 
of  Operator  Equations  (Lindberg;  Van  Rosendale,  Skeel) 

A  technical  report  with  this  title  was  written  by  Lindberg 

(1976).  The  abstract  follows: 

A  method  for  estimation  of  the  global  discretization  error 
of  solutions  of  operator  equations  is  presented.  Further 
an  algorithm  for  iterative  improvement  of  the  approximate 
solution  of  such  problems  is  given.  The  theoretical 
foundation  for  the  algorithms  are  given  as  a  number  of 
theorems.  Several  classes  of  operator  equations  are 
examined  and  numerical  results  for  both  the  error  estimation 
algorithm  and  the  algorithm  for  iterative  improvement  are 
given  for  some  classes  of  ordinary  and  partial  differential 
equations  and  integral  equations. 

The  table  of  contents  is  as  follows : 

1.  Introduction 

2.  General  Theory 

2.1.  Preliminaries 

2.2.  Basic  Theorems 

3.  Approximation  of  Linear  Functionals 

4.  Applications 

4.1.  Initial  Value  Problems  for  Ordinary  Differential  Equations 

4.2.  Two-Point  Boundary  Value  Problems  for  Ordinary  Differential 
Equations 

4.3.  Two-Dimensional  Elliptic  Boundary  Value  Problems 

4.3.1.  Problems  Nonlinear  in  u  Only 

4.3.2.  The  Minimal  Surface  Equation 

4.4.  Parabolic  Partial  Differential  Equations 

4.4.1.  The  Method  of  Lines  with  Euler's  Method 

4.4.2.  The  Method  of  Lines  with  the  Backward 
Euler  Method 

4.5.  Hyperbolic  Partial  Differential  Equations 

4.6.  Integral  Equations 

5.  Concluding  Remarks 


An  improved  and  shorter  version  of  this  report  (Lindberg  (1980))  is  to 
appear  in  BIT. 
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2.  A  Theoretical  Framework  for  Proving  Accuracy  Results  for  Deferred 
Corrections  (Skeel;  Ortman) 

A  manuscript  with  this  title  has  been  submitted  for  publication 

by  Skeel  (1980).  The  abstract  follows: 

General  techniques  are  described  for  proving  accuracy 
results  for  deferred  correction  solutions  to  differ¬ 
ential  equations.  These  techniques  apply  also  to 
computational  estimates  of  the  local  discretization 
error.  The  proofs  avoid  the  necessity  of  demonstrat¬ 
ing  the  existence  of  asymptotic  expansions  of  the 
global  error  in  powers  of  some  meshsize  parameter. 

The  outline  is  as  follows: 

1.  Introduction 

2.  Historical  Survey 

2.1.  Difference  Correction 

2.2.  Difference  Estimates  of  the  Local  Error 

2.3.  Generalized  Difference  Correction 

2.4.  Defect  Estimates  of  the  Local  Error 
2.5  "Cheap"  Global  Error  Estimates 

3.  Local  Error  Properties  of  Numerical  Methods 

4.  An  Error  Analysis  for  One  Deferred  Correction 

The  experiment  referred  to  in  the  penultimate  sentence  of 
Section  2.1  of  this  manuscript  was  performed  in  double  precision  on  the 
CYBER  175  for  a  decreasing  sequence  of  stepsizes  At  *  1,  1/2,...,  1/128 
for  the  problem  y'  =■>  y,  y(0)  =  1.  The  values  of  the  computed  solution  for 
t  *  1  were  printed  out,  and  they  seemed  to  be  converging  very  much  like 
(At)4. 


The  third  paragraph  of  this  manuscript  mentions  three  notable 
Implementations  of  the  deferred  correction  idea.  In  addition  there  is  the 
subroutine  SEPELI  developed  by  John  Adams  at  NCAR  (see  Swarztrauber  and 
Sweet  (1979)),  which  has  the  option  of  obtaining  fourth-order  accurate 
solutions  via  deferred  corrections  applied  to  a  fast  separable  elliptic  PDE 
solver  of  second  order  accuracy.  This  kind  of  situation  in  which  we  have 
a  very  efficient  low  order  method  is  ideal  for  the  application  of  deferred 


corrections. 


3.  The  Order  of  Accuracy  for  a  Deferred  Corrections  Algorithm  (Skeel) 

This  work  will  constitute  the  future  paper  mentioned  in  Skeel 
(1980,  section  4,  second  paragraph)  in  which  we  give  an  error  analysis 
for  a  sequence  of  iterations  for  the  algorithm  considered  by  Christiansen 
and  Russell  (1979) .  This  is  intended  to  be  a  realistic  example  of  the 
application  of  the  theoretical  framework  for  proving  accuracy  results, 
which  because  of  its  length  was  not  included  in  Skeel  (1980) . 

In  Christiansen  and  Russell  (1979)  a  careful  analysis  of  deferred 
corrections  that  does  not  involve  asymptotic  expansions  is  done  for  a 
realistic  algorithm  similar  to  the  implementation  of  Lentini  and  Pereyra 
(1977)  of  iterated  deferred  corrections  for  two-point  boundary  value 
problems.  Under  weak  assumptions  they  prove  that  each  iteration  increases 
the  order  by  one,  and  they  give  empirical  results  showing  that  the  first 
two  iterations  increase  the  order  by  two.  They  suggest  how  this  might  be 
proved  but  do  not  follow  through  because  "Such  a  proof  would  be  quite 
tedious."  In  this  paper  we  sketch  a  proof  of  this  fact,  which  we  believe 
is  less  tedious  than  that  of  Christiansen  and  Russell  (1979)  due  to  the 
way  in  which  we  break  down  the  proof  into  smaller  simply  stated  results. 

The  key  idea  is  the  use  of  judiciously  chosen  "discrete  Sobolev"  norms 
for  measuring  the  smoothness  of  the  global  error  and  of  the  local  error. 

In  particular  we  use  special  norms  like  those  of  Spijker  (1971)  and 
Stummel  (1975)  for  the  errors.  These  norms  are  chosen  so  that  they  admit 
both  upper  and  lower  bounds  on  the  norms  of  the  global  error  thus  making 
it  feasible  to  determine  the  exact  order  of  convergence.  However,  it 
must  be  acknowledged  that  the  algorithm  we  analyze  is  a  bit  simplified  and 
unrealistic  for  nonlinear  problems  in  that  we  assume  the  exact  solution  of 
nonlinear  equations. 
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3.1.  The  numerical  method.  Assume  that  the  differential  operator  F 
given  for  an  arbitrary  function  z  by 


F(z)  := 


b(z (0) ,  z(l)) 


yz' (x)  -  f (x,z(x)),  0  <  x  <  1, 
satisfies  the  assumptions  of  Christiansen  and  Russell  (1979)  so  that  in 
particular  the  operator  equation  F(y)  =  0  has  an  isolated  solution  y. 
Consider  a  family  of  meshes 

0  =  xo  <  xi  <••<  Xj  =  i 


with  hj  :=  x^  -  Xj_^  an<^  h  !=  1/J  such  that  h^  £  Ch  uniformly  for  all 
meshes  in  the  given  family.  (Presumably,  the  average  meshsize  h  can  be 
arbitrarily  close  to  zero,  for  otherwise  the  theorems  that  follow  have  no 
content.)  We  obtain  results  for  three  progressively  stronger  assumptions 
on  the  family  of  meshes: 

(i)  "no  assumption,"  meaning  that  there  is  no  assumption 
apart  from  that  already  stated, 

(ii)  "weak  assumption,"  meaning  that 
J-l 

I  |h  /h  -l|<ch 
j=l  J 

uniformly  for  all  meshes,  which  is  quite  realistic.  Skeel 

and  Jackson  (1979)  use  the  term  "variation-bounded"  to 

describe  such  a  family  of  meshes  in  connection  with  the 

stability  of  multistep  methods  for  initial  value  problems. 

(iii)  "strong  assumption,"  meaning  that 

max  |h  /h  -  l|  <  ch 
l<j<J-l  31  J 

uniformly  for  all  meshes,  which  is  not  so  realistic. 

The  term  "locally  uniform"  is  quite  descriptive  of  such 
a  family  of  meshes. 


W'.wV  •  **- 
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We  seek  a  numerical  solution  n 
theoretical  solution  on  the  mesh  y(x^), 
will  be  an  approximation  of  order  2k  to 


j 

j 


,  j  *  0(1)J,  approximating  the 
=  0(1) J.  Our  discretization 


g(y(xQ>,  y(xj))  =  0 


it. 

j  x 


j 

/  y'(x)  -  f(x,y(x))dx  =  0  ,  j  =  1(1)J  . 

j-l 


With  the  j-th  subinterval  [Xj_^,  x^  ]  we  associate  the  set  of  integers  I(j) 
consisting  of  those  2k  integers  from  0  through  J  which  are  nearest  j  -  1/2. 
We  approximate  the  term  f(x,y(x))  in  the  integrand  by  the  polynomial  of 
degree  2k-l  which  interpolates  f(x,y(x))  on  the  meshpoints  x^,  i  £  I(j). 
Thus  centered  interpolation  is  used  everywhere  except  near  the  boundary. 

The  resulting  equation  is 


T“-  (y(x  )  -  y(x  ))  Z  8.  .  f(x  ,y(x  ))  =  0(h2k) 
hj  j  J"1  i€l(j)  J*1  1  1 


where  g  is  a  function  of  the  relative  meshsizes. 
J  •  i 

Tlj  is  the  solution  of  the  discrete  problem 


The  numerical  solution 


g(nQ»  Hj)  =  0  , 


n .  ,)  -  I  8.  i  f(x.,  n.)  =  o  ,  j  =  1(1)J  , 
j_1  iel(j)  J,i  11 


which  is  a  system  of  J+l  nonlinear  equations  whose  Jacobian  has  bandwidth 

2k  times  the  dimensionality  of  the  original  ODE.  Thus  the  cost  of  solving 

2 

linear  systems  involving  the  Jacobian  is  proportional  to  k  .  This  can  be 
avoided  through  iterated  deferred  corrections. 

Let  X,  ■  (Cqi  Sj,...,  Cj)  be  an  arbitrary  element  of  the  discrete 
space  and  define  a  discrete  operator  <{>(£)  hy 

^k^O  *  8^0’  ’ 

■  h^  (cj "  Cj-i) "  i€j(J)  6j,i  f(v  V  » j  ■ 1(1)J  • 


For  k  *  1  this  is  the  trapezoid  rule 


‘MV  "  hj  tj-l)  "  2  f(xj*  V  “  2  f<xJ-l*  CjV  ’ 

which  is  the  most  economical  for  computation.  In  iterated  deferred 

corrections,  due  to  Fox  (1947),  we  write 

<J>  »:  4>  -  where  <t>:= 

mm  1 

thus  expressing  4>  as  the  sum  of  the  economical  trapezoid  rule  plus  the 
m 

correction  -  ^  .  Instead  of  solving  <f>(n)  -  (n)  “0,  the  procedure 

m  m 

is  to  solve 

•Kn1)  =  o 

and 

<Knk)  -  if>  (n1^1)  =  o  ,  k  =  2,  3 . 

m 

It  can  be  shown  that  the  high  order  corrections  are  unnecessary  for  small 
k  and  it  is  less  work  to  do  iterative  updating  deferred  correction : 

Hr\h  =  0  , 

4>(nk)  -  \(nk_1)  =  o  ,  k  -  2  (i)m  . 


3.2.  Stability.  We  are  interested  in  the  norm  of  the' error  ||f|k  -  Ay ||q 


where  Ay  =  (y(Xg),  y(x^),...,  y(Xj))  and  the  norm  is  defined  by 

IUIL*  max  |?  |.  We  can  reduce  this  to  the  local  level  because  <{>  is  stable 
V  0<j<J  3 


II?  -  cIIq  <  s0  !!♦«)  -  ♦<« ILq 


for  some  independent  of  the  mesh.  And  so 

llnk  -  Ay IIq  <  s0  ||<*>(nk)  -  4>(Ay)||*0  . 

The  norm  for  the  discrete  residual  space  E*^  has  not  yet  been  specified 
but  the  best  choice  is  the  Spijker  norm 
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,  J 

|!p)|  max  |p  +  E  hp| 

°  0<j<J  °  1“1 

where  p  Is  an  arbitrary  element  of  the  discrete  residual  space.  (This  norm 
also  yields  a  lower  bound  in  the  stability  definition.)  Stability  depends 
crucially  on  the  solution  being  an  isolated  solution  (Keller  (1976)). 

If  the  analysis  of  the  error  is  pursued,  it  turns  out  that 
accuracy  improvement  depends  on  the  smoothness  of  the  error.  For  this 
purpose  we  define  the  following  norms  for  the  discrete  solution  space: 
lid  :*=  max  (|?  |,  max  |tc,|)  ,  • 

l<j<J  J 

INI 2  :-max{|e0l,  ItcJ.  It^J) 

z  u  ±  2<j<J  2 


where 

^ j  (?j  “  Cj-l)/(xj  ”  Xj-1) 

and 

*-  (fcj  *  tCj_1)/(xj  -  Xj _2 )  . 

Define  the  following  norms  for  the  discrete  residual  space: 

||p||  max  {|pj,  max  |p  |}  , 

1  °  1<)<J  2 

||p||  max  { |  PQ  | ,  |p.|,  max  |fp  |}  '. 

2£j<J  J 


With  considerable  effort  one  can  establish  the  following  stability  results: 
(SJ  ||?  -  <  Sn  ||+(»  -  <K?)||*m,  m  -  0,  1,  2  . 

The  case  m  =  1  is  implied  by  Christiansen  and  Russell  (1979,  lemma  1)  for 
linear  problems. 


3.3.  Local  errors.  Returning  to  the  norm  of  the  error,  we  have  using  the 
more  general  norm 
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llnk  -  4>llm£  sm||*(nk)  -  *(Ay)||,o 

■  s»  H*t(''k"1)  -  • 

k  k _ 2. 

Thus  the  accuracy  of  n  depends  on  the  accuracy  of  )  as  an  estimate 

of  the  local  error  4>(Ay).  To  underscore  this  point,  we  write  the 
algorithm  as 

♦(n1)  -  o  , 

4»(nk)  -  ^(n1'-1),  k  =  2,  3 . 

the  idea  being  that  if  the  right  hand  side  were  exactly  the  local  error 
k 

then  ti  would  be  exactly  Ay.  The  error  in  the  local  error  estimate  can 

k-1 

be  split  into  the  propagated  error  arising  from  the  error  n  -  Ay  and 
the  local  error  arising  from  <j>^: 

*k(nk_1)  -  <|><Ay)  -  \(nk_1)  -  <|»k( Ay)  +  tpfe(Ay)  -  <KAy) 

-  [«k(nk_1)  -  \(Ay)]  -  4>k  (Ay)  . 

Each  of  these  terms  is  analyzed  separately;  in  this  subsection  we  examine 
-  <|>k(Ay) ,  which  is  the  local  error  of  the  formula  of  order  2k. 

With  "no  assumptions"  we  have  both 

(cj>  ll^CAy)!! *0  £  h2k  cj  ,  k  -  1,  2,... 

and 

(ck)  l!$k<Ay)l|*1  <  h2k  cj  ,  k  -  1,  2 . 

*2 

However  this  is  not  true  for  the  norm  for  k  ^  2  even  for  uniform  mesh 
families  because  the  use  of  uncentered  formulas  at  the  two  ends  causes 
abrupt  changes  in  the  local  error  there.  With  the  "strong  assumption"  one  can 


show  that 
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<c2)  l|4>(Ay)||*2  <  h2  c2  . 

Under  weaker  assumptions  this  Is  not  true  because  if  the  meshsize  does 
not  vary  slowly  then  neither  does  the  local  error. 


3.4.  Contractivity.  The  other  term  in  our  error  analysis  is 
|| -  ^(Ay)  11^.  If  there  is  to  be  an  increase  in  the  order  of 
accuracy,  the  operator  ^  needs  to  be  contractive  with  contractive  power 
at  least  0(h).  As  an  example. 


♦2<5)j 


i 

12 


(  h1+2hJ+1 

Vi+Wi 


t2f(xj,^j) 


V2Vi 

hj-i+hj+hj+i 


We  note  that  at  best 

-  o<\\z  -  ci^) 

so  that  contractivity  is  only  0(1).  However,  if  we  use  the  norm  for  the 

o 

discrete  solution  space  E~  involving  second  order  divided  differences,  then 

|*2 <Dj  -  *2(C)jl  =  0(h2  l|?  "  ^  * 

This  second  "bound"  is  better  than  the  first  if  £  is  smooth.  We  shall 
prove  contractivity  results  of  two  types: 

n*k«>  -  v5)ikihLLn!  -  cil 


and 

Note  that  more  contractivity  is  possible  if  we  are  willing  to  go  to  a 
stronger  norm. 

The  polynomial  interpolant  used  to  derive  and  hence  ^  can  be 
expressed  in  Newtonian  form  in  terms  of  divided  differences.  One  would 
find  that  for  1  £  j  <  J 
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j 

h  E  iMO, 

1-k  K 


j  i+k-1 
E  E 
i*k  S,-I-k+2 


ai,Jl  ^f(x£’  V 


j-1  i+k-1 

E  E 

i-k-1  £=i-k+2 


°i+l , £+1 


j+k-1 

E 

i-j-k-2 


J.* 


ff(x9>  C  )  -  h2 


2k-2 


ii  Vi+i 


tf(x,,  C.) 


Now  it  is  enough  to  show  that  n  ^f(Xj,  satisfies  a  Lipschitz  condition 
of  the  given  form.  For  the  first  Lipschitz  condition  this  follows  readily 
from  the  Lipschitz  continuity  of  f.  For  the  second  inequality  note  that 


tf<Y  +  V  -  *f(Y  sj> 

‘  VY  !3  +  eVd9'V 

•  1  tVYY*Vd<Kj + 1  V  Yi’c3-i+9Vi)d9’nj  ■ D 


Remark.  This  theorem  actually  holds  under  the  "weak  assumption." 
It  may  be  possible  to  further  weaken  the  assumption  on  meshes  for  the  first 
inequality  by  noting  that  only  a ^  -  2«i+1  £+1  +  ai+2  £+2  need  be  sma11 
rather  than 


“it  ”  °i+l,4+l  ’ 

THEOREM  2.  The  following  aontraativity  results  hold  for  the  E^1 
norm  with  "no  assumption" : 

a„>  ||*k(?)  -  *k(«||n  <  M*kl  ns  -  dij  . 

<>  ll*k«>  -  volln  i  h2li2  lit  -  dl2  - 

Proof.  It  is  sufficient  to  show  that  satisfies  the  given 

Lipschitz  conditions  and  this  we  do  as  in  the  proof  of  Theorem  1.  □ 

Remark.  It  may  be  possible  to  combine  a  contractivity  result  for 
^  with  a  stability  result  for  $  in  order  to  establish  a  stability  result 
for  <frk. 
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3.5.  Order  of  convergence.  We  separately  consider  the  results  for 
"strong  assumption,"  "weak  assumption,"  and  "no  assumption." 

With  the  "strong  assumption"  it  immediately  follows  from  the 
results  of  sections  2,  3,  and  4  that 

lln1  -  Ay  Ik  <  S2h2c2  , 


II H2  -  Ay  1^  <  S1(h2LJ2  Hn1  -  Ay^  +  h4cj)  , 

lln3  -  Ay Hq  <  S0(h2L^  ||n2  -  Ay +  h6cj)  , 
and  for  k  *  4 ,  5 ,  6 , . . . 

II nk  -  Ay ||q  <  S0(hL^0  ||nk_1  -  Aylk  +  h2kcj) 


so  that  the  order  progression  is  2,  4,  6,  7,  8,  9,....  The  result 

actually  proved  for  this  algorithm  by  Christiansen  and  Russell  is  that 

the  order  progression  is  2,  4,  5,  6,  7,  8»...  although  they  give 

empirical  evidence  for  the  stronger  result  and  outline  the  proof  for 

2 

this  result  under  the  stronger  assumption  that  h^+^/h^  =  1  +  0(h  ). 
With  the  "weak  assumption"  we  have 


lln1  -  Ay  1^  <  S^cJ  , 

lln2  -  Ay  Up  <  S0(hh201  lln1  -  Ay  Ik  +  h4c2)  , 
and  for  k  -  3,  4,  5,  6,... 


II nk  -  Ay H,  <  SqOiLqq  ||nk_1  -  Ay||0  +  h2kCk) 

so  that  the  order  progression  is  2,  4,  5,  6,  7,  8.... 
With  "no  assumption"  we  have 

lln1  -  Ay |k  <  S1h2c1  , 

and  for  k  -  2,  3,  4,  5,  6,... 

II  nk  -  Ay  Ik  <  S1(hLj1  II  nk_1  -  AyUx  +  h2kcj) 
so  that  the  order  progression  is  2,  4,  5,  6,  7,.... 
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4.  Ten  Ways  to  Estimate  Global  Error  (Skeel) 

A  manuscript  with  this  title  is  being  prepared  for  publication. 

An  outline  follows: 

0.  Introduction.  We  describe  and  assess  global  error 
estimation  techniques. 

1.  Deferred  correction. 

2.  Linearized  deferred  correction. 

3.  Differential  correction. 

4.  Linearized  differential  correction. 

5.  Defect  correction.  Conditions  are  established  sufficient 
for  the  validity  of  the  defect  correction  idea  of 

P.E.  Zadunaisky,  and  it  is  shown  that  this  approach  offers 
no  theoretical  advantage  over  deferred  correction. 

6.  Richardson  extrapolate  'n. 

7.  Error-gradient  estimation.  This  recent  if’ea  is  due  to 
Epstein  and  Hicks  (1979). 

8.  Using  two  different  tolerances. 

9.  Using  two  different  methods. 

10.  Using  a  method  with  a  specially  chosen  form  of  error. 

This  idea  of  Stetter  (1971)  has  a  weakness  which  has 
been  described  to  us  by  F.T.  Krogh. 
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5.  Nonsmooth  Situations  (Skeel) 

Ve  consider  the  use  of  a  single  deferred  correction  to  obtain  an 
improved  solution  and  hence  an  estimate  of  the  global  error.  Recall  that 
there  are  three  components  to  such  a  computation: 

(i)  a  numerical  solution  n  which  approximates  the  restriction 
to  a  mesh  Ay  of  the  theoretical  solution  to  the  operator 
equation  F(y)  =  0. 

(ii)  a  discretization  4>  of  F  which  is  computationally  attractive. 

(iii)  a  local  error  estimator  ip  for  <{>. 

A  corrected  solution  n  is  obtained  from  4>  (rj)  =  xp (rj)  •  The  theoretical 
justification  for  this  procedure  is  valid  only  under  certain  smoothness 
assumptions  on  the  original  solution  n  and  on  the  problem  F.  Numerical 
experiments  have  demonstrated  the  power  of  this  technique  when  the 
assumptions  are  satisfied,  but  when  they  are  not,  the  quality  of  the 
estimates  can  deteriorate,  in  particular,  the  order  of  accuracy  predicted 
by  the  theory  will  not  materialize.  This  section  discusses  in  general 
terms  what  might  be  done  in  nonsmooth  situations.  Application  to  parabolic 
PDEs  and  to  hyperbolic  PDEs  is  the  subject  of  sections  6  and  7. 

Of  special  interest  are  hyperbolic  PDEs  with  shocks,  and  for 
these  only  very  low  orders  of  accuracy  seem  possible  by  any  method  with 
the  possible  exception  of  very  complicated  shock  fitting  methods.  This 
Indicates  that  it  would  be  difficult  to  construct  asymptotically  correct 
estimates  for  the  local  error  near  discontinuities  in  the 
solution.  Hence  accuracy  considerations  for  error  estimation  techniques 
must  go  beyond  the  order  of  accuracy.  For  good  estimates  in  smooth 
regions  one  would  probably  desire  an  asymptotically  correct  estimate 
as  h  -*■  0.  Among  all  such  estimators  we  would  like  to  select  an  (inexpensive) 
estimator,  which  is  still  somewhat  accurate  when  the  solution  is  nonsmooth. 
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It  is  not  entirely  clear  how  to  measure  the  error  so  that 
|| n  -  Ay ||q  is  small  for  accurate  solutions,  since  a  slightly  displaced 
shock  gives  a  large  local  value  of  the  global  error  e  :=  n  -  Ay.  It 
may  be  better  to  use  a  root-mean-square  or  mean  norm  instead  of  the  max 
norm  for  ||  *  ||q*  Also,  it  might  help  to  define  Ay  as  average  values  on 
cells  rather  than  point  values,  as  in  the  "control  volume  approach" 
described  by  Roache (1975)  for  deriving  finite  difference  schemes.  For 
second  order  schemes  the  order  is  the  same  regardless  of  whether  we 
are  trying  to  compute  average  values  or  point  values,  but  it  does 
matter  for  higher  order  schemes. 

A  careful  examination  of  the  theoretical  basis  of  the  error 
estimation  technique  is  done  by  Skeel  (1980)  with  the  aim  of  identifying 
the  key  assumptions  necessary  for  the  success  of  this  technique.  In 
particular,  assumptions  involving  asymptotic  expansions  in  the  gridsize 
h  were  avoided.  It  was  found  that  the  discrepancy  ii  -  Ay  in  the  global 
error  estimate  satisfies  the  error  bound 

lln  -  Ay||0  <  S(hpK  ||e|^  +  ||i|>(Ay)  -  <J>(Ay)||*0) 
where  S  is  the  stability  constant  in 

lit  -  cIJq  <.  s  ||<fr<?)  “ 

and  Khp  is  the  contractivity  constant  in 

(*)  ||^P(t)  -  *KC)||*0  <  hPK  ||£  -  C|^  . 

Thus  there  are  four  factors  that  affect  the  accuracy  of  the  global  error 
estimate.  The  first  factor  is  the  norm  of  the  error  Uel^,  which  is 
defined  in  terms  of  the  first  q  divided  differences  of  e.  It  would  be 
desirable  to  modify  the  original  numerical  solution  r)  before  estimating 
the  local  error  so  as  to  reduce  the  norm  of  the  error.  This  may  be 
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possible  If  one  knows  something  about  the  behavior  of  £,  for  example  that 
e  alternates  in  sign.  The  second  factor  is  the  contractivity  h^K  of  the 
local  error  estimator  ip.  Very  often  ip  is  chosen  to  be  <p  ~  $  where  $  is 
a  more  accurate  discretization  of  F.  In  such  cases  one  should  choose  $ 
to  be  close  to  4>  in  some  sense.  The  third  factor  is  the  accuracy  of  the 
local  error  estimator  lJj(Ay)  -  <j>(Ay).  The  fourth  factor  is  the  stability 
constant  S  of  the  method  <}>..  In  the  remainder  of  this  section,  we  discuss 
each  of  these  in  more  detail. 


Accuracy  of  the  original  solution.  One  of  the  contributions 
to  the  inaccuracy  of  the  local  error  estimate  comes  from  the  original 
global  error  and  its  first  q  divided  differences.  The  local  error 
estimate  will  not  be  accurate  if  the  error  is  not  a  smooth  "function"  of 
the  independent  variables.  There  are  potentially  numerous  reasons  for 
roughness  in  the  global  error: 

(i)  roundoff  error.  This  would  normally  be  insignificant 
for  those  problems  of  interest. 

(ii)  iteration  error.  This  arises  from  the  use  of  implicit 
difference  schemes  and  can  be  minimized  by  doing 
enough  iterations. 

(iii)  the  use  of  more  than  one  difference  scheme  to  approximate 
the  differential  equation.  This  situation  occurs  when 
two-level  difference  schemes  are  used  to  calculate 
starting  values  for  three-level  schemes;  it  also  occurs  in 
hyperbolic  equations  when  implicit  schemes  are  used  to 
calculate  numerical  boundary  values  for  explicit  schemes, 
(lv)  coarseness  of  the  grid.  Certain  problems  such  as  stiff 
differential  systems  can  be  accurately  solved  even  though 
the  grid  is  too  coarse  for  the  asymptotic  theory  to  hold. 
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(v)  irregularly  spaced  gridpoints. 

(vi)  nonsmoothness  in  the  problem  iteslf  including  the  boundary. 

If  something  is  known  about  the  behavior  of  the  error,  it  is 
possible  by  smoothing  the  numerical  solution  T)  to  reduce  the  error  and 
especially  the  divided  differences  of  the  error.  For  example,  in  the 
case  of  shock  calculations,  higher  order  schemes  introduce  systematic 
oscillatory  errors  and  random  choice  methods  introduce  random  errors  into 
the  solution.  In  our  operator  notation  this  means  the  decomposition 

•  X  into  a  smoothing  operator  x  and  a  simpler  local  error  estimator 
$.  Experimentation  on  a  simple  test  problem  with  1/4,  1/2,  1/4  smoothing 
in  space  resulted  in  error  reduction  near  the  shock  but  not  away  from  the 
shock.  It  should  be  possible  to  refine  such  techniques  by  the  use  of  a 
switch  triggered,  for  example,  by  sign  changes  in  second  differences  of  the 
numerical  solution.  A  related  approach  is  used  by  Lindberg  (1976)  for 
the  leap-frog  method.  There  the  local  error  estimate  for  a  gridpoint  was 
constructed  from  computed  values  only  at  every  other  gridpoint.  Either 
of  these  approaches  is  applicable  to  stiff  ordinary  differential  equations 
solved  by  the  trapezoidal  rule.  If  less  is  known  about  the  error  behavior, 
one  can  treat  the  roughness  as  noise  and  construct  approximations,  such  as 
least-squares  approximations,  which  filter  out  the  noise. 

Contractivity  of  the  local  error  estimator.  One  basis  for 
evaluating  potential  local  error  estimators  is  inequality  (*).  One 
determines  the  weakest  differentiability  assumptions  and  the  smallest 
numerical  constant  h^K.  This  can  also  be  done  for  reduced  Integer  values 
of  p  and/or  q.  Under  weaker  smoothness  assumptions  it  is  quite  acceptable 
that  the  contractivity  of  ip  be  some  fraction  of  unity  rather  than  0(h). 
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It  should  also  be  fruitful  to  study  the  contractivlty  for 

simple  test  problems,  in  particular,  the  constant  coefficient  linear 

problem.  This  can  be  used  to  determine  the  contractivlty  of  ip  or  better 

still  the  contractivlty  of  ip  •  <p  \  since  the  Frechet  derivative  of  the 

later  operator  maps  the  original  local  errors  into  inaccuracies  in  the 

local  error  estimate.  Following  is  an  example  of  this  type  of  analysis 

in  a  situation  where  asymptotic  analysis  in  powers  of  h  is  inappropriate: 

Example.  Consider  a  stiff  system  of  ODEs  y'  =  f(t,  y)  with 

initial  conditions.  Let  $  be  the  finite  difference  operator  for  the 

backward  Euler  method.  Define  the  global  errors  £  :=  n  -  y(t  )  and 

n  n  n 

the  local  errors 

50  :=  °»  6n  :=  "  y(tn-l))/h  +  f(tn’  y(tn})  ’ 

Consider  the  local  error  estimators 

An)  7  (f(t  ,  n  )  -  f(t  .,  n  ,))/h 

n  x  n  n  n— l  n-x 


and 


/(n)„  :=  ^  f^^(t  ,  n)  . 


n  n 


Then  it  can  be  shown  that  for  the  test  problem  y1  =  Xy  +  g(t).  Re  X  £  0, 

/(Ay  +  e)  -  /(Ay)  =  ^  l  (1  -  hX)n-j+1  (6.  -  6  )  -:  p  , 

n  n  i  j=1  j  j-1  n 

which  is  small  compared  to  6^,  since 

|Pnl  <in"1/2|«il  +1  |(6  -  6  )/h|  . 

2<j<n  J  J 


However , 


/(Ay  +  e)  -  /(Ay)n  «p  -fi, 
n  n  n  z  n 


which  can  be  very  large  compared  to  6 
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Accuracy  of  the  local  error  estimator.  A  potential  estimator 
4>  can  also  be  evaluated  by  determining  the  weakest  differentiability 
assumptions  and  the  smallest  numerical  constant  for  a  bound  on 
||^(Ay)  -  <j>(Ay)||A0.  This  can  also  be  done  for  lower  orders  of  accuracy. 

Stability  of  the  method.  For  global  error  estimation  it  may 
be  worthwhile  to  consider  methods  <p  for  the  "improved"  solution  which 
are  different  from  the  original  method  $.  For  example,  for  conservation 
laws  one  might  choose  a  good  monotone  method  for  <p  since  it  may  respond 
in  a  more  stable  fashion  to  the  perturbations  <Kn)  that  are  introduced 
into  the  second  integration. 
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6.  Parabolic  Problems  (Skeel;  Van  Rosendale) 

A  couple  of  nonsmooth  problems  of  this  type  were  tested  by 

Lindberg  (1976) :  the  heat  equation 

u  -  u 
t  xx 

with  nonsmooth  initial  conditions  and  Burgers'  equation 

u  +  uu  =  Vu 
t  X  XX 

with  a  small  value  for  V.  This  second  problem  develops  a  sharp  front 

with  large  spatial  derivatives.  The  local  error  estimators  of  Lindberg 
E  E  E 

have  the  form  ^  :=  <p  -  <p  where  $  is  a  higher  order  discretization  of 

the  operator  based  on  polynomial  interpolation.  We  attempted  to  improve 
his  results  by  using  a  higher  order  discretization  $  having  a  more 
compact  computational  molecule  more  like  that  of  $.  This  would  hopefully 
yield  a  more  contractive  $:=<}>-$. 

Statistics  were  generated  for  both  local  and  global  error 
estimates  at  each  meshpoint  in  space-time.  These  included  not  only  the 
correct  error  and  the  inaccuracy  of  the  error  estimate  but  also  the 
two  separate  contributions  to  the  inaccuracies:  the  part  inherited  from 
the  error  of  the  original  numerical  solution  and  the  part  introduced  by 
the  local  error  estimator.  These  two  contributions  were  evaluated  by 
splitting  the  inaccuracy  of  the  local  error  estimate  <Kn)  ~  4>(Ay)  into 
tp(n)  “  'KAy)  plus  -  <p(Ay).  Perhaps  it  would  have  been  better  to  split 
the  inaccuracy  of  the  global  error  estimate  according  to 

n  -  Ay  -  (n  -  n)  +  (ri  -  Ay) 

where  fi  solves 


<Kn)  -  <KAy)  . 

Thus  we  could  have  separately  evaluated  the  centractivity  of  ip  given  by 
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|| ?i  -  Ti||/||n  -  Ay ||  and  the  relative  error  of  ip,  which  is  in  some  sense 
measured  by  ||n  -  Ay||/||n  -  Ay  }|.  This  information  might  be  meaningful 
even  if  ip  is  only  to  be  used  for  local  error  estimation. 

Results  for  our  compact  local  error  estimators  are  better 
than  for  Lindberg's  estimators  in  regions  of  smoothness  but  are  as  bad 
at  discontinuities. 


6.1.  Heat  equation.  The  differential  equation 

u  —  u  =0 
t  XX 

was  solved  with  zero  boundary  conditions  and  with  both  smooth  initial 
conditions 

u(x,  0)  =  sin  x  ,  0  £  x  £  tt  , 

and  nonsmooth  initial  conditions 

u(x,  0)  =  min  {x,  ir-x}  ,  0  £  x  £  tt  . 

Analytic  solutions  are  given  by  Lindberg  (1976). 

The  original  numerical  solution  U  was  obtained  by  FTCS  differencing 


i  6  U^+l/2  _  _L  $2  U?  =  0 


k  t  “j 


2  VI  w 

*  i 


where  u”  ~  u(jh,  nk).  This  is  second  order  in  the  spatial  gridsize 

2 

h  *  1/J  provided  that  the  stability  condition  k  £  h  /2  is  satisfied. 
This  scheme  was  also  used  as  the  basic  method  (f>  for  obtaining  the 
corrected  solution. 

E  EE 

Lindberg  (1976)  considers  ip  :m  <p  -  <p  where  $  is  based  on 
quadratic  interpolation  in  time  and  quartic  interpolation  in  space 

i  Mt«t  -  0 

with  modifications  for  j  ■  1,  j  -  J-l  and  n  ■  1.  We  also  consider  the 
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operator  $  given  by 

1  n  .  1  ,2. .  nn+l/2  _1_  .2  n+1/2  n 

k  (1  +  12  6x)6t  Uj  "76xytUj  "  °  * 

because  it  is  the  method  fourth  order  in  h  which  differs  the  least  from 

E  E 

<p.  The  local  error  estimator  ip  :=$-</>  is  given  for  arbitrary  V  by 

k  1  ,2  n  h2  1  j.4  „n 

2^'t  Vj  -T2^6xVJ 

and  :=  $  -  $  is  given  by 

,k  h2.  1  x  .2  n+1/2 

We  are  interested  in  how  contractive  these  operators  are.  Since 
they  are  both  linear,  it  is  sufficient  to  consider  ||l(»^(V)||  and  ||<KV)|| 
where  the  norm  is  the  max  norm.  We  obtain  the  bounds 


II*K(V)II  <  +  |v|t 


and 


11*00  II  < 


h  -  6k 


2 

3n  k 


min  l\-  IVI^,  k|v|t) 


“2  2  n 

where  |  v| ^  is  the  maximum  of  all  |h  Vj  I  an<*  1^1 1  the  niaximum  of 

,-l 


all  |k  6t  v"|.  One  can  show  that 


(1) 


bound  for  ip  < 
E  — ■ 

bound  for  ^ 


h  -  6k 


3h  +  2k 


s\ 


assuming  that  k  £  n  /2.  With  even  weaker  norms  we  get 


ll*E<v)  ||  <  ||v|| 

3h  k 


11*00  111 


2h  -  12k 


3h2k 


llvll  , 


and 
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and  again  (1)  holds  for  this  last  pair  of  bounds.  This  suggests  that  $  is 

more  contractive  than  <j>E.  It  is  worth  noting  that  $  becomes  identical  to 
2 

<J>ifk  =  h  /6  with  the  result  that  the  error  estimates  would  be  zero! 
Nonetheless,  zero  would  be  a  good  approximation  to  the  error  in  an 
absolute  error  sense . 

The  table  that  follows  compares  for  the  smooth  problem  the 

error  estimates  at  x^  for  a  uniform  grid  0  =  Xq  <  x^  <. . .<  x^  =  ir 

2  2  E  “ 

with  k  =  (5/2ir  )h  .  The  estimated  global  errors  are  U  -  U  and  U  -  U 

E  EE**  - 

where  U  solves  <j>(U  )  =  (U)  and  U  solves  4»(U)  =  iKU)i  The  errors  in 

g 

the  local  error  estimates  ip  (U)  and  ^(U)  can  be  split  into  two  parts: 

<//E(U)  -  <J>(Au)  =  [ipE(Au  +  e)  -  i{jE(Au)  ]  -  4>E(Au)  , 

ip(U)  -  <KAu)  =  [i|j(Au  +  e)  -  ip (Au) ]  -  $(Au)  . 

The  first  part  is  due  to  the  error  e  :=  U  -  Au  in  the  original  solution  U, 

and  the  second  part  is  the  local  error  introduced  by  the  more  accurate 
E  £  — 

operators  <f>  =  <P  -  ip  and  <p  *  <p  -  ip. 


time 

level 

global 

error 

error  in  global 
error  estimate 

local 

error 

error  in  local  error  estimate 
*  part  due  to  e  +  new  error 

for  \J)E 

for  ty 

for  i|)E 

for 

1 

2 

4 

8 

-157 

0 

0 

-5083 

-10 

135  -145 

-10 

7  -17 

-304 

-15 

-1 

-4929 

-464 

-771  307 

-9 

7  -16 

-572 

-40 

-1 

-4634 

-434 

-723  289 

-7 

8  -15 

-1011 

-77 

-1 

-4096 

-381 

-637  255  • 

-3 

—12 _ si L 

Note:  all  values  are  10^  times  actual  values. 

The  next  table  displays  the  same  information  for  the  nonsmooth 


problem: 
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time 

level 

global 

error 

error 

in  global 
estimate 

local 

error 

error  in  local  error  estimate 
*  part  due  to  e  +  new  error 

error 

,  .E 
for  ip 

for  ip 

for  ipE 

for  <p 

1 

-177 

139 

139 

-5741 

4500 

12355  -7855 

4500 

9294  -4793 

2 

-147 

54 

110 

-745 

-1705 

251  -1956 

133 

161  -28 

4 

-114 

51 

81 

-308 

-295 

-595  300 

41 

48  -7 

8 

-:84 

47 

58 

-116 

-52 

-102  50 

10 

11  -2 

A 

Note:  all  values  are  10  times  actual  values. 

Hence  for  both  problems  ip  was  a  much  more  accurate  local  error 
£ 

estimator  than  ip  although  the  big  local  error  near  the  point  of 
discontinuity  was  still  underestimated  by  a  factor  of  4.6.  The  better 
global  error  estimator  was  ij;  due  to  a  fortuitous  cancellation  of  errors  in 
the  local  error  estimates  resulting  from  the  use  of  a  special  local  error 
estimate  for  the  first  time  level.  Apparently  the  two  parts  of  the  error 

in  the  local  error  estimate  tend  to  cancel  especially  in  the  case  of  the 

smooth  problem.  It  is  not  clear  whether  or  not  this  desirable  situation 
can  be  caused  to  happen  more  generally.  The  foregoing  experiments  were 
performed  also  for  h  =  tt/19  and  h  =  it/ 39 .  The  qualitative  nature  of  the 

results  seemed  to  depend  on  the  time  level  n  rather  than  the  time  t. 

The  previously  tabulated  calculations  were  also  performed  for 
the  uniform  grid  0  =  x^  <  x^  <...<  x^q  *  it  with  k  »  (5/2ir^)h^.  Error 
estimates  are  given  at  x^,  which  lies  right  on  the  discontinuity: 
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time 

level 

global 

error  in  global 
error  estimate 

local 

error  in  local  error  estimate 
*  part  due  to  E  +  new  error 

error 

for  i{iE 

for  \l> 

error 

for 

for  ip 

1 

193 

-330 

-330 

7703 

-13218 

-31213  17995 

-13218 

-26542  13324 

2 

146 

-519 

-213 

-1235 

-12344 

-29350  17006 

-90 

-129  40 

4 

106 

-203 

-148 

-420 

-735 

-1262  527 

-1 

8  -9 

8 

76 

-120 

-104 

-142 

-95 

-161  66 

0 

3  -2 

4 

Note:  all  values  are  10  times  actual  values. 

We  note  that  the  global  error  estimates  are  really  bad  due  to  the  very 
poor  local  error  estimate  at  the  initial  discontinuity.  Otherwise  ip  was 

E 

much  more  accurate  than  ip  . 

Clearly  a  nonuniform  mesh  would  be  appropriate  for  the  nonsmooth 
problem,  and  so  the  grid  points 

*3  ■  f  t(4(j  -  i>2  +  1>(j  -  f>  +  1J 
were  used  with  J  =  10  and  the  time  increment 


k  -  4  (min  x  -  x  n)2 
j  3  3 


The  gridpoint  density  is  four  times  as  great  at  the  center  as  at  the  ends. 
The  method  <j>  is  given  by 


«  ufl/2  -  jL  tfV  .  o 

{t  1  Sx2  3 


where 


6  ..n+1/2  _  1  ,„n+l  ..n 

6t  Uj  k  (Uj  "  Uj 


and 
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1  "H  ^7  <“"«  -  *k»  • 

which  is  generally  only  first  order  in  space.  The  higher  order  method 
$  is  given  by 

.  1_  .  .  62  .  1  ^  .  v  A  ^  A  „n+1/2  -i?-  i^1/2  -  n 

(!  +  12  hjhj+1  ^  +  3  (hj+1  -  h^)  6x  )  6t  Uj 


6x 


where 


and 


uj+1/2  *  i  (uj  +  uj+1) 


^ ^  (v  <u”+1  •  ^  (u5  -  ”J-1»  • 

which  is  generally  only  third  order  in  space.  Hence  the  local  error 
estimator  ip  :=  $  -  $  is  given  by 


,k  _l_  Ail  Vth-1/2  I  r*  nil  vn+l/2 

(2  “  12  hjhj+l)  6t  6x2  Vj  ”3  (hj+l  V  St  6x  ;j 

Numerical  results  are  given  below  for  x^  =  ir/2. 


time 

level 

global 

error 

error  in  global  error  estimate 
for 

local 

error 

2 

120 

-133 

-434 

8 

82 

-59 

-135 

20 

80 

-36 

-35 

30 

83 

-29 

-19 

Note:  all  values  are  10  times  actual  values. 

There  is  an  improvement  compared  to  the  uniform  mesh  due  to  the  smaller 
contribution  to  the  global  error  from  the  local  error  at  the  discontinuity. 
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6.2.  Burgers*  equation.  Burgers*  equation  u^  +  uu^  “  Vu^  w*t*1  Initial 
condition  u(0,  x)  ■  sin  x,  0  <  x  <  ir,  and  zero  boundary  conditions  was 
solved  with  an  implicit  scheme 


♦ 


n 

n-1 


Xj-1  Xj  Xj+1 


which  is  first  order  in  k  and  second  order  in  h.  Two  difference  operators 
were  considered  as  local  error  estimators: 


j-2 


o 

-O - 0 - 0 - 0 

6 

xj-l  xj  xj+l  xj+2 


t 

n 


n-1 


n-1 


with  modifications  for  j  ■  1,  j  *  J-l,  and  n  *  1;  and 


$ 


Xj-1  XJ  Vi 


E 

The  difference  operator  4>  ,  which  is  second  order  in  k  and  fourth  order  in 

h,  was  used  by  Lindberg;  and  <j>,  which  is  second  order  in  both  k  and  h,  was 

chosen  because  it  differed  the  least  from  <f>.  The  following  tables  compare 

the  two  local  error  estimates  with  h  =  tt/20  and  k  -  1/10  for  V  *  1  and 

for  V  -  1/16.  In  this  problem  a  front  develops  and  then  dissipates.  The 

sharpness  of  the  front  increases  as  v  +  0.  From  the  table  it  appears  that 
E  — 

<P  is  better  than  $  except  at  the  first  time  level  where  unsymmetric 
formulas  had  to  be  used.  The  tables  give  the  maximum  absolute  value  of 


each  quantity. 
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V  -  1 


time 

level 

global 

error 

error 

error 

In  global 
estimate 

local 

error 

error  in  local 
=  part  due  to 

error  estimate 
e  +  new  error 

for 

E 

for  iJj 

for  i{j 

E 

for  ip 

1 

107 

35 

8 

1483 

704 

510 

199 

134 

113 

118 

5 

188 

32 

20 

399 

74 

47 

29 

139 

95 

43 

10 

195 

20 

20 

211 

30 

25 

7 

45 

32 

14 

20 

143 

11 

4 

73 

4 

2 

2. 

1 

2 

2 

40 

39 

1 

1 

10 

0 

1 

0 

1 

1 

0 

A 

Note:  all  values  are  10  times  actual  values. 


V  **  1/16 


time 

level 

global 

error 

error  in  global 
error  estimate 

local 

error 

error  in  local 
*  part  due  to 

error  estimate 
e  +  new  error 

for  ipE 

for  ip 

-  ,  E 

for  ip 

_ 

for  ip 

1 

45 

12 

2 

478 

183 

140 

•  42 

27 

49 

41 

5 

184 

45 

11 

398 

179 

112 

67 

68 

82 

62 

10 

262 

53 

205 

962 

462 

132 

448 

824 

1611 

806 

20 

1839 

1039 

2478 

5179 

1508 

3620 

2112 

5529 

10946 

5417 

40 

1054 

340 

931 

1104 

401 

179 

580 

1063 

2188 

1125 

80 

213 

15 

129 

72 

27 

26 

2 

66 

136 

70 

Note:  all  values  are  10  times  actual  values. 
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7.  Hyperbolic  Problems  (Skeel;  Ortman,  Van  Rosendale) 

Guided  by  the  general  theory  of  section  2  and  the  Ideas  of 
section  5,  various  numerical  schemes  were  constructed  and  tested. 

The  estimation  of  global  discretization  error  by  deferred 
corrections  requires  two  parallel  integrations  of  the  problem,  the  second 
integration  being  corrected  by  subtracting  out  local  error  estimates 
obtained  from  the  first  integration.  Thus  there  are  three  major  components 
of  any  numerical  test:  the  initial  value  problem,  the  integrator,  and 
the  local  error  estimator. 

Recently  Stetter  (1978)  has  observed  that  for  global  error 
estimation,  the  second  integration  could  be  performed  by  the  cheapest 
available  (stable  and  consistent)  method.  However,  there  are  economies 
associated  with  re-using  the  same  method,  and  so  it  is  not  clear  in 
general  which  approach  is  more  efficient.  In  all  our  tests  the  two 
integrations  are  performed  with  the  same  numerical  method. 

Also,  we  note  that  Hackbusch  (1977)  and  authors  cited  therein 
have  applied  Richardson  extrapolation  to  methods  for  hyperbolic  systems 
of  first  order. 

Three  ways  of  constructing  the  local  error  estimator  ip  are 
outlined  in  Skeel  (1980,  section  4).  We  were  attracted  to  the  Fox-Stetter- 
Lindberg  idea  of  using  ip  :=  <p  -  <p  where  d>  is  the  basic  integrator  and  $  is 
a  more  accurate  discretization.  This  permits  the  use  in  ip  of  specialized 
techniques  for  nonlinear  hyperbolic  problems  such  as  artificial  viscosity, 
upwind  differencing,  artificial  compression,  antidiffusion,  and  hybridization. 
It  is  worth  noting  that  as  a  solution  operator,  $  need  not  be  stable  in 
order  for  it  to  be  useful  for  local  error  estimation.  In  addition,  the 
computational  cost  of  "inverting"  $  is  irrelevant  since  only  $(n)  need 
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be  computed.  Hence,  certain  types  of  implicitness  are  computationally 
inexpensive. 

Three  initial  value  problems  with  known  solutions  were  tested: 

(i)  u^  +  u^  =  0  with  a  periodic  square  wave  solution,  (ii)  the  inviscid 

2 

Burgers  equation  u  +  (u  /2)  =  0  with  initial  conditions  chosen  so  that 

t  x 

a  shock  develops,  and  (iii)  the  Riemann  problem  for  the  system  of  three 
equations  for  one-dimensional  Eulerian  gas  dynamics  as  in  Sod  (1977). 

The  first  method  <p  considered  was  the  two-step  Lax-Wendroff 
scheme.  A  variable-mesh  scheme  was  desired,  but  it  was  shown  that  there 
is  no  second  order  extension  of  Lax-Wendroff  to  variable  meshes.  A 
reasonably  compact  two-level  fourth  order  error  estimator  for  Lax-Wendroff 
was  sought  but  only  third  order  estimators  could  be  found  that  were  also 
economical.  Finally,  smoothing  schemes  for  Lax-Wendroff  solutions  were 
tested  with  partial  success.  The  aim  was  to  make  the  numerical  solution 
more  amenable  to  global  error  estimation. 

Also  tested  was  Lax-Wendroff  with  artificial  viscosity  as  in 
Sod  (1977). 

Since  most  methods  have  a  fractional  order  of  convergence  in 

the  norm  when  applied  to  problems  with  shocks,  it  was  thought  worthwhile 

to  test  first  order  difference  schemes.  Deferred  correction  applied  to 

FTCS  (forward  time  centered  space)  differencing  with  a  Lax-Wendroff  local 

error  estimator  yielded  a  solution  which  was  better  than  that  obtained  by 

Lax-Wendroff  alone.  When  Lax's  method  was  teamed  up  with  Lax-Wendroff,  the 

deferred  correction  solution  was  better  than  the  original  Lax  solution 

but  not  as  good  as  Lax-Wendroff.  The  problems  tested  were  ut  +  u^  «  0 

and  u  +  uu  ■  0. 
t  X 
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A  number  of  techniques  for  hyperbolic  systems  of  conservation 
laws  seem  to  be  motivated  by  considerations  appropriate  to  a  single 
conservation  law.  These  involve  nonsmooth  functions  of  differences  such 
as  absolute  value,  sign,  and  maximum.  Systematic  extension  of  these 
schemes  to  systems  +  ^(u)x  =  .0  is  possible:  one  way  is  to  perform  a 
local  transformation  which  diagonalizes  _f ' (u)  at  some  local  value  of  u 
making  the  system  only  weakly  coupled  locally,  then  to  apply  the  scalar 
scheme,  and  finally  to  reverse  the  transformation.  These  transformations 
can  be  done  implicitly.  For  example,  for  upwind  differencing  this  might 
amount  to 

,i  <^+1  -  ij) 

*  j  sgn  (f '  (u" ) )  (f(Uj+1)  -  2f(uJ)  +  f<uJ_j)) 

where  sgn(f' (11))  is  the  matrix  obtained  by  diagonalizing  _f'  (u) ,  applying 
the  sign  function  individually  to  each  eigenvalue,  and  reversing  the 
similarity  transformation.  This  idea  is  used  by  Lax  and  Wendroff  (1960) 
for  artificial  viscosity,  and  they  show  how  to  do  it  given  the  eigenvalues, 
but  not  the  eigenvectors,  of  _f ' (u) .  Usually  the  eigenvalues  of  _f ' (u) 
are  known,  and  in  any  case,  there  is  always  the  possibility  of  approximating 
sgn^f'  0j))  without  reducing  the  order  of  accuracy.  We  tested  the  foregoing 
upwind  scheme  on  the  Euler ian  gas  dynamics  equations,  and  it  performed  much 
better  than  the  upwind  differencing  tested  by  Sod  (1977).  This  idea  may 
be  useful  for  artificial  compression,  antidiffusion,  artificial  dissipation, 
and  SHASTA. 

Deferred  corrections  can  also  be  used  as  a  means  of  obtaining 
more  accurate  solutions.  For  hyperbolic  equations  it  could  be  viewed  as 
a  type  of  hybrid  method  which  combines  the  desirable  error  propagation 
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properties  of  a  less  accurate  method  with  the  smaller  local  errors  of 
a  more  accurate  method.  For  example,  the  Glimm-Chorin  random  choice 
method  does  an  excellent  job  of  resolving  shocks,  but  it  is  at  best 
first  order  accurate.  Hence  the  solution  obtained  may  benefit  from 
smoothing  followed  by  a  deferred  correction  with  some  second  order 
method.  Another  example  is  upwind  differencing,  which,  because  it  is 
monotone,  always  converges  (Harten,  Hyman,  and  Lax  (1976))  to  the 
physically  relevant  solution  of  a  scalar  PDE.  Since  it  is  first  order, 
it  could  also  benefit  from  a  deferred  correction.  It  can  be  shown 
that  if  one  combines  two  conservative  methods  in  this  way,  then  the 
resulting  deferred  correction  method  is  conservative. 


7.1.  Test  problems.  Three  hyperbolic  conservation  laws  were  used  with 
initial  conditions  chosen  so  that  the  analytical  solution  is  known. 

The  first  problem  is 


_a_ 

3t 


u 


0 


with  initial  condition 

u(x,  0)  = 

and  periodic  boundary  condition  u(l,  t)  =  u(0,  t) .  The  analytic  solution 
is  a  periodic  square  wave  travelling  to  the  right  at  unit  velocity. 

The  general  nonlinear  scalar  conservation  law  has  the  form 

■h  u  +  £<u)  ■  0  • 


For  initial  condition  u(x,  0)  *  g(x)  this  can  be  solved  via  the  hodograph 
transformation  to  yield  the  implicit  analytical  solution 
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u  =  g  (x  -  f '  (u)  t) 

for  u  =  u(x,  t).  The  solution  to  the  differential  equation  ceases  to 
exist  at  points  where  u  is  about  to  become  multiple-valued.  But  a  weak 
solution  exists  having  a  shock  discontinuity  originating  at  such  points. 

The  speed  S  =  S(t)  of  such  a  shock  satisfies  the  Rankine-Hugoniot  condition 

f(u_)  -  f(uT) 


UR  "  UL 

where  uT  =  u(S(t)-,  t)  and  uD  *  u(S(t)+,  t).  However  not  all  discontinuities 
L  K 

satisfying  this  condition  are  realizable  from  continuous  initial  conditions. 

Imagine  the  jump  discontinuity  u  -  u  being  split  into  u  -  u  and 

K  L  L 

Uj^  -  u  where  u  is  strictly  between  u^  and  u^.  In  order  for  these  two 
shocks  not  to  separate,  we  must  have 

f(u)  -  f(uL)  f(uR)  ”  f(u) 


u  - 


“L 


UR  -  U 


This  inequality,  which  must  hol<  for  all  u  strictly  between  u^  and  uR, 

is  known  as  the  entropy  condition. 

1  2 

The  choice  f(u)  =  ^  u  corresponds  to  the  inviscid  Burgers 
equation,  whose  solution  satisfies 

u  ■  g(x  -  ut)  . 

The  initial  condition  u(x,  0)  =  g(x)  can  be  chosen  to  give  interesting 
analytic  solutions.  For  example,  one  can  arrange  to  have  two  shocks 
develcp  spontaneously  and  then  later  merge  into  a  single  shock.  We  used 
the  initial  condition  u(x,  0)  =  -arctan  x  and  imposed  boundary  conditions 
at  x  =  +  1  obtained  from  the  analytic  solution  which  satisfies 
u  ■  -arctan (x  -  utj  and 


0  <  u  <  v/2  for  x  <  0  , 
u  *  0  for  x  =  0  , 


-ir/2  <  u  <  0  for  x  >  0  . 
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For  t  >  1  there  is  a  stationary  shock  discontinuity  at  x  =  0. 

For  a  system  of  M  conservation  laws 

^  ICu)  =  0 

one  might  be  able  to  obtain  a  solution  to  the  Riemann  problem,  which 
has  the  initial  conditions 


u(x,  0)  = 


X  <  X 
X  >  X 


0  ’ 
o  • 


For  this  problem  we  hypothesize  a  piecewise  constant  splution  with 

values  Uq,  u^,...,  u^  separated  by  M  discontinuities  at  x  =  x^  +  S^t, 

m  =  1(1)M.  The  weak  form  of  the  PDE  reduces  to  the  Rankine-Hugoniot 

conditions  at  the  discontinuities: 

S  (u  -  u  .)  =  f(u  )  -  f(u  -),  m  *  1(1)M  . 
m  — m  — m— 1  —  “in  —  — tn— i 

2  2 
This  represents  M  nonlinear  algebraic  equations  for  the  M  unknowns 

— 1’  1*  ant*  Sl’  Sm"  ^are  must  be  taken  to  fix  up  the 

solution  at  those  discontinuities  where  the  entropy  condition  is  not 

satisfied. 


The  equations  for  the  iulerian  specification  of  one  dimensional 


gas  dynamics  are 


where  p  is  density,  u  is  velocity,  p  is  pressure,  and  e  is  internal 
energy  per  unit  mass.  We  assume  a  perfect  gas  for  which  the  equation 
of  state  is 


P  -  (Y  -  l)pe  . 

The  specific  heat  ratio  y  is  chosen  to  be  7/5,  which  is  its  value  for 
diatomic  gases  such  as  air.  It  is  sometimes  convenient  to  use 
conservation  variables 


where  m  is  momentum  and  e  is  energy  per  unit  volume.  If  the  fourth 
variable  is  eliminated  by  the  equation  of  state,  then 


l(u)  =f  (Y  -  De  +  \  (3  -  Y)®2/P 
^yme/p  -  j  (y  -  l)m3/p2 


The  analytic  solution  of  the  Riemann  problem  for  these  equations  reduces 
to  a  single  nonlinear  equation,  which  can  be  solved  by  the  Godunov 
iteration.  At  discontinuities  where  the  entropy  condition  is  violated, 
one  must  use  a  rarefaction  wave  solution  instead.  The  initial  conditions 
used  are  those  of  Sod  (1977)  : 

at  t  =  0  we  have  p=l,  u=0,  p=l  for  0  £  x  <  y  and 

p  »  1/8,  u  -  0,  p  •  1/10  for  y  <  x  <  1  . 


7.2.  Two-step  Lax-Wendroff  method.  For  the  problem  u^  +  f(u)^  =  0  this 
method  consists  of  a  step  of  Lax's  method 

*72  -  \  >w>  +  5  sx  ‘OW:- 0 

followed  by  a  leapfrog  step 


~  6  u"+1/2  +  ~  5  f (u"+1/2)  -  0  , 

k  t  j  h  x  j  * 


best  illustrated  schematically  by 


«  leap 

frog . 

3 - 

[  Lax  j 

] - 

[  Lax  | 

For  an  arbitrary  function  v  the  Lax-Wendroff  finite  difference  operator  <(• 
yields 
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n+1  n+1/2  k  n+1/2  -  - 

4>(Av)j  *  (vt  +  f(v)x)j  -  j  (f'(v)(vt  +  f(v)x))x  +  0(h  +  k  )  . 

With  v  ■  u  we  get  the  local  error 

<|>(Au)J+1  -  h2  ¥  +  0(h3) 

where 

¥  (cu  -  u  )  -  4  (f '  (u) (cu  -  u  )) 

24  tt  xx  t  8  tt  xx  x 

and 

c  :=  k/h  . 

The  global  error  satisfies 

-  uj  -  h2e"  +  0(h3) 

where  e  =  e(x,  t)  solves 

et  +  (F’(u)e)x  -  -  ¥  . 

Only  the  expansion  for  (f>(Av)  is  needed  for  developing  deferred  correction 
error  estimators. 

We  managed  to  obtain  fairly  close  agreement  with  the  numerical 
results  of  Sod  (1977)  for  the  Lax-Wendroff  solution  of  the  Riemann 
problem.  For  the  time  increment  he  chose 

k  =  0.9  max(|u|  +  c)h  . 

For  the  analytic  solution  the  maximum  is  attained  in  Region  4  of  Figure  3 
of  Sod.  From  c  =  /yp/ P  and  Table  II  of  Sod  we  get  k  =  0.411h.  The  value 
of  t  in  Figures  4  through  15  is  not  given,  but  it  is  clear  that  the  shock 
position  is  very  close  to  .75.  From  equation  (15)  and  Table  II  the  shock 
speed  can  be  calculated  and  we  get  t/k  ■  35.  For  the  Lax-Wendroff  method 
an  artificial  viscosity  term  was  added.  In  our  test  we  also  ran  the 
problem  without  artificial  viscosity,  and  in  spite  of  greater  oscillations, 
the  average  (^  norm)  error  at  t  -  35k  was  reduced  by  10%  for  the  density. 


24%  for  the  pressure,  43%  for  the  velocity,  and  16%  for  the  internal 
energy  with  similar  reductions  at  earlier  time  levels.  (Note:  the 
column  labelled  "e"  in  Table  II  should  be  labelled  "e"  and  the  figures 
labelled  "ENERGY"  should  be  labelled  "INTERNAL  ENERGY.")  This  may 
reflect  a  shortcoming  of  the  SL ^  norm, which  does  not  sufficiently  penalize 
oscillations.  It  has  been  suggested  that  the  & ^  norm  may  be  more 
appropriate  because  convergence  results  have  been  obtained  for  this  norm. 

For  the  nonlinear  scalar  test  problem  a  variable  mesh  is 
certainly  appropriate,  and  for  this  reason  we  sought  a  second  order 
accurate  variable-mesh  generalization  of  the  Lax-Wendroff  method.  In 
order  to  avoid  problems  near  boundaries  and  to  limit  computational 
costs,  it  was  required  that  a  linear  combination  of  u”, 

f(U?),  U“  and  f(U?,,)  and  that  be  a  linear  combination  of 

j  j+i  j+l  3 

u",  f(u"),  f(uj+l)*  and  f(Uj+l/2)*  The  method 


,.n 


U?  +  (x 


-  xj  u" 


j+1/2  x  +1  -  x  lVAj+l  ~j+l/2'  uj  VAj+l/2  uj+l 

-  \  (P(0j+1)  -  F(U"»] 


of1  -  U?  -  - 

3  j  x. 


(f<u”S/2>  -  F<uj:i/2» 


j+l/2  -  j  —1/2 

satisfies  these  restrictions  provided  that  the  whole-numbered  meshpoints 
satisfy 

X j  "  2  (xj-l/2  +  Xj+l/2)  ; 

otherwise,  there  is  no  second  order  formula.  Given  the  whole-numbered 


meshpoints,  it  is  always  possible  to  solve  for  the  half-numbered  meshpoints 
but  it  may  not  be  possible  to  force  them  to  lie  between  the  appropriate 
whole-numbered  meshpoints.  It  would  be  better  to  define  the  half-numbered 
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meshpoints  to  be  midway  between  the  whole-numbered  meshpoints  and  to 
define  new  whole-numbered  meshpoints  to  be  midway  between  the  half- 
numbered  meshpoints. 

It  is  well  known  that  without  enough  artificial  viscosity 
Lax-Wendroff  produces  oscillations  next  to  a  shock.  For  a  bounded 
stationary  numerical  solution  to  ut  +  f(u)^  =  0  where  f(u)  is  strictly 
convex  or  concave  one  can  show  that  to  a  first  approximation  the  error 
decreases  geometrically  with  alternating  sign  as  one  moves  away  from 
the  stationary  shock.  In  section  5  we  indicated  that  errors  can  be 
better  estimated  if  they  are  smooth.  A  numerical  experiment  was 
performed  for  the  Lax-Wendroff  method  without  artificial  viscosity 
applied  to  the  inviscid  Burgers  equation  with  h  =  0.1  and  k  =  0.05. 

The  table  that  follows  contains  values  for 

error  in  u"  , 


t  1  ,.n  ■  1  T _n  *  1  „n 

error  In  j  ?  Uj  +  j  Dj+J 


error 


ln  4  '  j-1/2  +  Uj+l/2  +  Uj-l/2  +  j+1/2^  ’ 


The  -jj.  smoothing  is  exact  for  linear  polynomials  and  removes 
alternating  errors  of  constant  magnitude.  The  second  smoothing  is 
suggested  by  the  simple  analysis  mentioned  previously.  If  we  look  at 
the  errors  in  the  worst  case,  it  is  clear  that  the  second  smoothing  is 
better  than  the  first  smoothing  which  in  turn  is  better  than  not 
smoothing  at  all.  Similar  results  were  obtained  for  h  ■  0.05  and  k  ■  0.025. 
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134 

• 

320 

• 

303 

• 

.0228 

.0742 

.0701 


.0008 

.0995 

.0927 


039 

““  • 

217 

• 

207 

• 

-.0067 

.0197 

.0195 


-.0174 

.0052 

.0078 


3 
2 

.0012 


-.0026 

.0011 

.0013 


-.0003 

.0009 

.0010 


-.0002 

.0008 

.0008 


-.2006 

.0763 

.0604 


i 

-.3692 

.0759 

“  • 

IS 

07 

-.0001 

1 

1.30 

.0454 

-.0541 

• 

24 

.0005 

.0166 

-.0083 

—  • 

£ 

05 

.0005 

-.0024 

-.0002 

.0010 

.0007 

.0011 

.0008 

-.0024 

-.0002 

.0011 

.  0007 

.0010 

.0007 

-.0025 

-.0002 

.0006 

.0006 

.0007 

.0006 

In  'the  remainder  of  this  subsection  we  consider  local  error 
estimators  f  «  ^  ^  for  the  Lax-Wendroff  method  4>.  We  begin  by 

restricting  ourselves  to  operators  $  involving  only  two  time-levels  n 
and  n+1.  Freely  available  are  values  of  f  at  u”,  ant*  for 

any  j.  The  order  of  $  cannot  exceed  the  order  of  i  in  time, which 
can  be  determined  by  sending  h  ■+•  0.  Thus  we  can  begin  our  search  for 
$  by  considering  discrete-time  continuous-space  operators,  which  are 
little  more  than  ODE  methods  for  IVPs.  The  Lax-Wendroff  method  becomes 

Un+1  .  u"  +  k  F(Un  +  F(Un)) 

where  F  denotes  the  operator  -f(»)x  and  Un  *  Un(x).  With  the  values  of 
F  at  Un,  U°  +  ^  F(Un),  and  Un+^  it  is  not  possible  to  construct  a 
third  order  implicit  Runge-Kutta  method  <j>.  If  we  permit  one  additional 
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^+1/2  "  Ux  Uj+l/2  "  2h  6x  f(Uj+l/2)  ’ 

Uj+1/2  =  (1  -  |  62)uJ  -  ^  6x  f(u“+1/2)  , 

TTn+ 1  » i  1  T7n  ^  »  r  - /TTn+ X  /  2  \ 

Uj+l/2  ~  Mx(1  "  8  6x)Uj+l/2  "  h  5x[f(Uj+l/2) 

+  lSx  ^+1/2' 1  • 


“T  - -  H  'f  £<i^+1)  + 1  <•*, + 1  - 1 5x>  <1/2> 

♦Iv'-K'  • 

We  compared  this  method  to  Lax-Wendroff  for  the  gas  dynamics  problem 
without  artificial  viscosity.  The  average  error  for  the  fourth  order 
method  at  t  =  35k  was  greater  by  8%  for  the  density,  20%  for  the  pressure, 
64%  for  the  velocity,  and  14%  for  the  internal  energy.  Also  it  produced 
a  54%  overshoot  in  the  velocity  compared  to  33%  for  Lax-Wendroff. 

If  one  is  willing  to  use  a  three-level  operator  <j>,  then  one  can 
avoid  additional  evaluations  of  f  by  using 


t  “tV1 + i  + h  ux6x(i + i  st>  f(ui+l)  •  °  • 

The  prospects  of  obtaining  good  error  estimates  with  these 
operators  $  did  not  seem  very  good,  and  so  we  turned  our  attention  to 
estimating  errors  of  first  order  methods. 


7.3.  Forward-time  centered-space  differencing.  This  method  is  known  to 
be  unstable  for  the  linear  model  problem,  but  it  was  accidentally  tested 
due  to  erroneous  programming  of  Lax's  method.  The  linear  test  problem 
was  solved  with  h  *  1/10  and  k  =  1/60,  and  a  deferred  correction  solution 
was  computed  using  Lax-Wendroff.  Below  we  have  the 
error  in  the  FTCS  solution 
error  estimate  from  deferred  correction 


8 


for  all  the  spatial  meshpolnts  at  t  ■  30k: 


30 

-.25 

-.69 

• 

o 

o 

1.06 

.00 

.13 

-.87 

.19 

.12 

31 

-.11 

-.46 

00 

o 

• 

1 

.47 

-.08 

.44 

-.47 

.09 

-.08 

Out  of  interest  we  also  computed  the  Lax-Wendroff  solution,  and  below  we 
have  the 

error  in  deferred  correction  solution 
error  in  Lax-Wendroff  solution 


• 

o 

o 

-.14 

-.22 

• 

o 

00 

.59 

00 

o 

• 

-.31 

o 

• 

1 

.11 

.21 

.12 

-.16 

-.34 

.07 

.74 

-.05 

-.13 

-.60 

.20 

.15 

In  sum,  deferred  correction  performed  well  in  this  instance. 

7.4.  Lax's  method.  This  is  a  highly  dissipative  first  order  scheme 
given  by 

“f  ’  Ha  *  “I+1>  *  I  "A  • 

For  the  linear  test  problem  with  h  =  1/10  and  k  =  1/60  we  tried  for  $ 
both  Lax-Wendroff  and  Lax-Wendroff  with  diffusive  and  antidiffusive 
fluxes  as  described  by  Sod  (1977).  In  both  cases  the  results  were 
disastrous  probably  due  to  the  fact  that  the  numerical. solution  separates 
into  two  uncoupled  solutions,  one  with  j+n  even  and  the  other  with  j+n 
odd.  Lindberg  (1976)  encountered  the  same  problem  with  the  leapfrog 
scheme.  Out  of  interest  we  also  computed  solutions  directly  with  each 
of  the  $  schemes  and  the  average  error  was  17Z  less  with  diffusive  and 
antidiffusive  fluxes. 

We  repeated  the  experiment  with  leapfrog  for  $  and 
with  k  ■  1/50.  The  deferred  correction  solution  was  nearly  the  same 
as  the  Lax  solution,  and  thus  was  useless  for  global  error  estimation. 
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In  order  to  avoid  separated  solutions  we  considered  instead 
a  two-step  Lax  scheme  in  which  we  first  compute  Uj^/2  using  the  Lax 
method  and  then  compute  from  these  two  values  again  using  the  Lax 

method.  With  k  =  1/60  and  $  =  Lax-Wendroff  we  computed  a  deferred 
correction  solution.  Below  we  have  the 
error  in  the  Lax  solution 
error  estimate  from  deferred  correction 
for  all  the  meshpoints  at  t  =  20k: 


.23 

.26 

.30 

.35 

-.63 

-.63 

.34 

.30 

.25 

.23 

.14 

.07 

-.02 

-.12 

-.15 

-.14 

-.07 

.02 

.12 

.15 

These  are  not  very  good  results.  The  same  combination  for  <j>  and  (J>  was 
also  tested  on  Burgers'  equation  with  h  =  1/10  and  k  =  1/50.  The  results 
were  much  better.  Below  we  give  the 

time  level 

average  error  in  Lax  solution 

average  error  estimate  from  deferred  correction 


1 

3 

7 

15 

30 

0014 

.0041 

.0099 

.0233 

.0572 

0013 

.0039 

.0090 

.0194 

.0417 

We  also  compared  the  deferred  correction  solution  to  the  Lax-Wendroff 
solution  and  found  that  by  the  30th  time  level  the  average  error  of  the 
latter  was  smaller  by  a  factor  of  14  but  that  spurious  oscillations  were 
less  pronounced  for  the  former.  The  deferred  correction  error  was 
smaller  than  the  Lax  error  by  a  factor  of  4. 


7.5.  Upwind  differencing.  For  Burgers'  equation  we  tested  the 
conservative  upwind  scheme 
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1  ,  n+1  n*  ,  1  ,  n  n  s  n 

k  (uj  "  Uj5  +  h  (gj+l/2  "  ®j-l/2  0 


where  (omitting  n) 

gj+l/2  =  2  f^Uj+l^  +  2  f*Uj*  "  2  s8n^f'('2  uj+l  +  2  Uj^f*Uj+l^  “  f(uj^  * 

If  f'(u)  is  positive  for  u  close  to  u^,  this  has  the  effect  of  defining 

u?+1  in  terms  of  u?  ,  and  u?,  which  is  desirable  because  the  characteristic 
j  J-l  3 

passing  through  (Xj ,  tn+^)  would  pass  between  (x^  tn)  and  (x^ ,  tn).  If 
f '  (u)  is  negative,  then  uj+^  given  in  terms  of  u”  and  uj+^*  With 
h  *  1/10,  k  =  1/50,  and  $  =  Lax-Wendroff  we  computed  a  deferred 
correction  solution.  Below  we  give  the 
time  level 


average  error  in  upwind  solution 

average  error  estimate  from  deferred  correction 


1 

3 

7 

15 

30 

35 

0003 

.0010 

.0025 

.0054 

.0106 

.0112 

0004 

.0011 

.0025 

.0053 

.0089 

.0092 

Out  of  interest  we  also  computed  the  Lax-Wendroff  solution,  and  below 
we  have  the 

time  level 

average  error  in  deferred  correction  solution 
average  error  in  Lax-Wendroff  solution 


1 

3 

7 

15 

30 

35 

00002 

.00007 

.00017 

.00051 

.00238 

.00335 

00002 

.00007 

.00017 

.00046 

.00164 

.00246 

For  the  gas  dynamics  equation  Roache  (1975,  p.  237)  defines 


upwind  differencing  by 
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nn+1  .  nn  .  k  rrn  r>n  \  k  /cn  e11  \ 

h  "j  -  s«"  h  <%  -  -  2h  <^1+1  -  Vl>  > 


;j+l  -j- 


2  T  T  ti 

where  G  =  (pu,  pu  ,  u(e+p))  ,  _S  =  (0,  p,  0)  and  sgn  =  sgn(u^).  (The 

equations  given  by  Sod  (1977,  p.  30)  appear  to  be  erroneous.)  The 

numerical  results  obtained  by  Sod  were  not  good.  There  was  a  nonphysical 

shock  where  there  should  have  been  an  expansion  wave.  Our  results  were 

somewhat  similar  except  that  there  were  also  oscillations  present.  The 

stability  condition  CT  £  1  given  by  Sod  differs  from  the  very  stringent 

condition  given  by  Roache.  This  latter  condition  is  violated  in  the 

experiments. 


Because  of  the  failure  of  this  upwind  differencing  scheme  for 
the  system  of  three  hyperbolic  equations,  we  tried  to  extend  the  scalar 
scheme  used  for  Burgers'  equation  to  a  system  of  conservation  laws.  A 
systematic  way  of  doing  this  was  discussed  in  the  introduction  to  this 


section.  For  upwind  differencing  it  involves  using  sgn(f' (u))  where  if 
(“)  =  Q  diag(X^,  X2,...,  Am)Q  \  then  sgn(f'  (u))  =  Q  diag(sgn(A1) , 
sgn(A2)»*«-»  sgn(AM))Q  ,  which  is  independent  of  Q.  For  the  gas 
dynamics  equations  it  is  most  convenient  to  use  the  variables  p,  u,  and 
c  *  /y (y-l)e,  and  after  substitution  into  (u)  we  have 


■  Q  diag(u,  u+c,  u-c)Q 


where 
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I 

I 


1 

u  +  c 


u  u  ,  .  c 

T  T  +  uc  +  =F 


l 

u  -  c 

2  2 
u  ,  c 

-r-  -  UC  +  — ; 


Y-l  "  2 

u 

2 

u  uc 

u  c 

4  2 (y-l) 

2  2 (Y-l) 

2 

u  .  uc 

u  c 

^4  +  2  (Y-l) 

‘  2  "  2 (Y-l) 

The  following  may  be  useful  in  organizing  the  computation: 


u  “  Q  *  t  I/2 


£(u)  -  f ' (u)u  . 

Remark.  The  nonconservative  upwind  differencing  given  by 

Roache  would  be  in  line  with  the  above  approach  if  we  had  used  the 

Y_  2  1  3 

splitting  _f  =  G  +  j>  where  £  =  (Pu,  pu  ,  P“)  as  l°ng  as  lul  <  c. 

The  matrix  upwind  differencing  scheme  described  in  the 
paragraph  before  the  remark  was  tested  on  the  Riemann  problem,  and  it 
produced  a  very  nice  solution.  With  $  =  Lax-Wendroff  (without 
artificial  viscosity)  we  computed  a  deferred  correction  solution. 
Compared  to  the  original  matrix  upwind  solution,  the  deferred  correction 
solution  had  average  errors  at  t  *  35k  which  were  worse  by  more  than  a 
factor  of  2  for  all  variables.  In  addition  a  spurious  rarefaction  shock 
and  compression  shock  appeared  between  the  rarefaction  expansion  and 
the  contact  discontinuity. 


I 


.1 

.  I 

I 

u 


vi 
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A  compact  second  order  but  implicit  scheme  is  the  box  scheme 
1  «  ITn+l/2  ,  1  r  , /TIn+l/2N  n 

k  V*  "j+1/2  +  h  4  A  £(Dj+l/2)  '  °  • 

Whether  this  should  be  used  to  determine  U1?^  from  or  vice  versa 

J+l  J 

depends  on  the  direction  of  the  characteristics.  A  conservative  upwind 
box  scheme  is  given  by 

„n+l  TTn  k  f  n+1/2 

U.-  =  U  -  —  o  g. 

j  j  h  x 

where 

n+1/2  h  .  »  ,Tn+l/2  ,  c/.Itn+l/2. 

®j+l/2  "  4k  4t4x  Uj+l/2  +  f(V/2> 

-i.gn 

+  "t4*  f<"jS/2>J  • 

If  we  were  to  use  this  as  a  solution  method,  one  could  apply  this 
iteratively  beginning  with  Uj+^  ®  for  all  j .  Doing  one  iteration  is 

equivalent  to  our  first  order  upwind  scheme,  and  doing  two  iterations 
yields  second  order  accuracy.  This  scheme  extends  to  a  system  as  we 
have  described.  We  tested  this  scheme  as  a  solution  method  for  the 
Riemann  problem  and  found  that  with  two  iterations  per  time  level  the 
average  error  in  the  solution  variables  was  comparable  to  that  of  the 
first  order  upwind  differencing.  With  three  iterations  per  time  level 
the  error  was  decidedly  worse. 

A  deferred  correction  solution  was  computed  from  the  upwind 
solution  of  the  Riemann  problem  using  $  ■  upwind  box  with  2  iterations. 
(We  discovered  that  <j>(U)  would  remain  unchanged  if  we  had  instead  used 
$  *  implicit  upwind  box.)  Below  we  have  given  for  three  different  time 
levels  the 
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average  errors  in  matrix  upwind  solution 
average  error  estimates  from  deferred  correction 


time  level 

P 

P 

u 

G 

5 

.0052 

.0024 

.0059 

.0032 

.0107 

.0082 

.0175 

.0127 

15 

.0088 

.0049 

.0082 

.0052 

.0152 

.0113 

.0294 

.0193 

35 

.0125 

.0075 

.0108 

.0071 

.0200 

.0139 

.0440 

.0269 

We  also  computed  the  solution  in  two  other  ways,  and  below  we  have  the 
average  errors  at  various  time  levels  for 

deferred  correction  with  $  =  upwind  box  with  2  iterations 
deferred  correction  with  $  =  upwind  box  with  3  iterations 
Lax-Wendroff  (without  artificial  viscosity) 


time  level 

P 

P 

u 

e 

.0043 

.0044 

.0087 

.0118 

5 

.0040 

.0040 

.0084 

.0108 

.0064 

.0069 

.0150 

.0181 

.0056 

.0052 

.0105 

.0192 

15 

.0051 

.0046 

.0093 

.0174 

.0079 

.0071 

.0156 

.0242 

.0073 

.0062 

.0148 

.0314 

35 

.0070 

.0059 

.0143 

.0299 

.0084 

.0064 

.0127 

.0279 
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II.  GAUSSIAN  ELIMINATION  AND  NUMERICAL  INSTABILITY 

The  solution  of  linear  systems  of  equations  is  basic  to  the 
numerical  solution  of  many  problems,  and  yet  this  subproblem  has  not 
been  treated  in  an  entirely  satisfactory  way.  Stewart  (1973)  states  that 
"In  spite  of  intense  theoretical  investigation,  there  is  no  satisfactory 
algorithm  for  scaling  a  general  matrix."  A  theoretical  solution  to  the 
scaling  problem  was  discovered  by  Skeel  (1979).  This  was  one  of  the 
results  of  a  study  of  the  implications  of  a  stability  concept  which  is 
more  appropriate  for  sparse  systems.  This  investigation  was  stimulated 
by  a  report  of  Gear  (1975). 
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1.  Scaling  for  Numerical  Stability  in  Gaussian  Elimination  (Skeel;  Ortman) 
A  paper  with  this  title  was  authored  by  Skeel  (1979).  The 
abstract  follows: 

Roundoff  error  in  the  solution  of  linear  algebraic 
systems  is  studied  using  a  more  realistic  notion 
of  what  it  means  to  perturb  a  problem,  namely,  that 
each  datum  is  subject  to  a  relatively  small  change. 

This  is  particularly  appropriate  for  sparse  linear 
systems.  The  condition  number  is  determined  for 
this  approach,  the  effect  of  scaling  on  the  stability 
of  Gaussian  elimination  is  studied,  and  it  is 
discovered  that  the  proper  way  to  scale  a  system 
depends  on  the  right-hand  side.  However,  if  only  the 
norm  of  the  error  is  of  concern,  then  there  is  a  good 
way  to  scale  that  does  not  depend  on  the  right-hand  side. 

The  table  of  contents  is  as  follows: 

1.  Introduction 

2.  Condition  of  Linear  Systems 

3.  Stability  of  Algorithms  for  Linear  Systems 

4.  Gaussian  Elimination  with  Column  Pivoting 

4.1.  Scaling  for  Numerical  Stability 

4.2.  Scaling  for  Accuracy 

5.  Gaussian  Elimination  with  Row  Pivoting 

5.1.  Scaling  for  Numerical  Stability 

5.2.  Scaling  for  Accuracy 

6.  Practical  Implications 

Appendix  A.  Error  Bounds  for  Column  Pivoting 
Appendix  B.  Error  Bounds  for  Row  Pivoting 

Measurements  of  certain  quantities  introduced  in  this  paper  were 
made  for  the  UNPACK  SGEFA  test  problems.  Unfortunately  for  our  purposes 
these  problems  were  not  representative  of  those  likely  to  occur  in 
practice,  but  rather  there  were  a  number  of  unusual  problems  designed 
to  test  the  logic  of  SGEFA.  For  all  but  one  of  the  problems  the  backward 
error 


max 


b  -  Ax] 

ww 


was  less  than  the  unit  roundoff  error  u  even  though  the  ill  scaling  ratio 
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2.  Iterative  Improvement  Implies  Numerical  Stability  for  Gaussian 
Elimination  (Skeel) 

A  paper  with  this  title  was  authored  by  Skeel  (1980).  The 
abstract  follows: 

Because  of  scaling  problems,  Gaussian  elimination  with 
pivoting  is  not  always  as  accurate  as  one  might  reasonably 
expect.  It  is  shown  that  even  a  single  iteration  of 
iterative  refinement  in  single  precision  is  enough  to 
make  Gaussian  elimination  stable  in  a  very  strong  sense. 
Also,  it  is  shown  that  without  iterative  refinement  row 
pivoting  is  inferior  to  column  pivoting  in  situations 
where  the  norm  of  the  residual  is  important. 

The  table  of  contents  is 

1.  Introduction 

2.  Numerical  Stability 

3.  Gaussian  Elimination  with  Column  Pivoting 

4.  Error  Bounds 

5.  Backward  Error  Bounds 

6.  Gaussian  Elimination  with  Row  Pivoting 


61 


3.  Effect  of  Equilibration  on  Residual  Size  for  Partial  Pivoting  (Skeel) 
A  paper  with  this  title  was  authored  by  Skeel  (198x).  The 
abstract  follows: 

It  is  shown  that  column  pivoting  with  row  equilibration 
satisfies  the  same  type  of  error  bound  as  does  row 
pivoting  without  scaling  and  that  row  pivoting  with 
column  equilibration  satisfies  the  same  type  of  bound  as 
does  column  pivoting  without  scaling.  An  interesting 
consequence  for  row  pivoting  is  that  column  equilibration 
is  sufficient  to  ensure  that  the  norm  of  the  residual  is 
reasonably  small. 

The  table  of  contents  is 


1.  Introduction 

2.  Main  Results 

3.  Proofs  of  Results 
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III.  NUMERICAL  SOLUTION  OF  ORDINARY  DIFFERENTIAL  EQUATIONS 

Systems  of  ODEs  occur  naturally  In  many  applications  and  also 
arise  from  the  spatial  discretization  of  time-dependent  partial 
differential  equations.  The  most  effective  numerical  methods  for  solving 
initial  value  problems  in  ordinary  differential  equations  have  been 
multistep  methods.  Work  in  this  area  is  reported  here. 
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1.  Equivalent  Forms  of  Multistep  Formulas  (Skeel) 

The  Adams-Bashforth-Moulton  and  the  backward  differentiation 

formulas  are  very  popular.  Improved  stability  properties  are  possible 

with  other  formulas,  but  they  are  not  used  for  several  reasons.  One  of 

the  difficulties  concerns  storage  requirements.  It  is  shown  in  Skeel  (1979) 

how  to  formulate  the  corrector  formula  and  how  to  select  a  predictor  so 

that  any  multistep  method  can  be  implemented  as  cheaply  as  the  two 

popular  ones.  The  abstract  of  this  paper  follows: 

For  uniform  meshes  it  is  shown  that  any  linear  k-step 
formula  can  be  formulated  so  that  only  k  values  need  to  be 
saved  between  steps.  By  saving  additional  m  values  it  is 
possible  to  construct  local  polynomial  approximations  of 
degree  k  +  m  -  1,  which  can  be  used  as  predictor  formulas. 
Different  polynomial  bases  lead  to  different  equivalent 
forms  of  multistep  formulas.  In  particular,  local  monomial 
bases  yield  Nordsieck  formulas.  An  explicit  one-to-one 
correspondence  is  established  between  Nordsieck  formulas 
and  k-step  formulas  of  order  at  least  k,  and  a  strong 
equivalence  result  is  proved  for  all  but  certain 
pathological  cases.  Equivalence  is  also  shown  for  P(EC)* 
formulas  but  not  for  P(EC)*E  formulas. 

The  table  of  contents  is  as  follows: 

1.  Introduction 

2.  Minimum  Storage  for  Multistep  Formulas 

3.  Construction  of  Linear  Nordsieck  Formulas 

4.  The  Correspondence  Between  Multistep  and  Nordsieck  Formulas 

5.  Equivalence  of  Linear  Nordsieck  Formulas  to  Linear 
Multistep  Formulas 

6.  Equivalence  of  Predictor-Corrector  Formulas 

7.  Applications 
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2.  The  Stability  of  Variable-Step  Nordsleck  Methods  (Skeel  and  Jackson) 
Work  on  the  stability  of  interpolatory  step  changing  has  been 
prepared  for  publication  by  Skeel  and  Jackson  (198x).  The  abstract 
follows : 


Conditions  are  given  for  the  stability  of  the  Nordsieck 
formulation  of  Adams  and  backward  differentiation  methods. 

It  is  shown  that  if  the  stepsize  selection  function  is 
variation-bounded,  then  these  methods  are  stable.  It  is 
proven  that  for  each  method  there  exist  constants  a,  b, 
c  satisfying  a  £  b  <  1  <  c  such  that  if  the  stepsize  ratio 
h  +1/h  satisfies  0  <  h  ./h  <  a  or  if  b  <  h  /h  <  c 
tRen  tRe  method  is  stab?e7  ?or  k- value  methoRs  wi?h  k  £  8, 
tables  are  given  for  the  values  a,  b,  and  c. .  If  the 
stepsize  is  kept  constant  for  one  or  more  steps,  then 
methods  are  more  stable  in  the  sense  that  the  intervals 
defined  by  a,  b,  and  c  are  larger.  Tables  and  graphs  are 
given  which  show  how  the  stability  of  a  method  changes  as 
the  stepsize  is  changed  less  often. 

The  table  of  contents  is  as  follows: 

1.  Introduction 

2.  General  Stability  Considerations  and  Theoretical  Framework 

3.  Stability  of  Adams-Bashforth-Moulton  Methods 

4.  Stability  of  Backward-Differentiation  Methods 

5.  Conclusion 
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3.  Blended  Linear  Multistep  Methods  (Skeel;  Dahlquist-Downs,  Vu) 

At  present  the  most  popular  codes  for  stiff  systems  of  ODEs 
are  based  on  the  backward  differentiation  formulas.  One  promising 
alternative  is  to  use  the  "blended"  multistep  formulas  (Skeel  and  Kong 
(1977))  which  attempt  to  combine  the  best  features  of  the  Adams  and  the 
backward  differentiation  formulas.  By  definition  every  autonomous 
stiff  system  must  have  at  least  one  "component"  which  is  nonstiff  and 
the  Adams  formulas  are  very  good  for  such  equations.  There  are 
theoretical  reasons  for  believing  that  the  blended  methods  are  both 
effective  and  computationally  inexpensive.  Limited  empirical  evidence 
suggests  that  the  blended  formulas  may  be  as  good  as  the  backward 
differentiation  formulas  for  stiff  problems,  better  for  nonstiff 
problems,  and  much  better  for  stiff  oscillatory  problems.  We  have  been 
working  on  implementing  these  methods  by  making  modifications  to  state- 
of-the-art  computer  codes  for  stiff  ODEs. 

In  1978  from  the  University  of  Toronto  we  acquired  STIFF  DETEST, 
which  is  a  testing  program  for  stiff  ODE  integrators.  We  attempted  to 
execute  the  program  on  an  INTERDATA  minicomputer,  but  the  object  code 
was  too  long,  and  so  we  wrote  our  own  testing  program.'  An  improved 
version  of  STIFF  DETEST  was  acquired  in  1979,  and  we  have  successfully 
executed  it  in  five  and  a  half  hours  on  a  departmental  PRIME  computer. 

It  seems  that  the  execution  time  is  much  too  great,  and  we  are  looking 
into  this  problem. 

The  blended  formulas  had  been  Implemented  in  a  structured 

version  of  new  DIFSUB,  which  has  been  distributed  to  several  researchers. 

We  compared  the  efficiency  of  this  to  that  of  GEAR,  Rev.  3  of  Hindmarsh 

(1974)  on  several  test  problems.  Testing  was  performed  as  described  by 

SI-  1  and  Kong  (1977)  except  that  we  used  the  CYBER  175  and  h  . 

min 

se*.  to  4utf,  where  u  is  the  unit  roundoff  error  2”^. 
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e 

max 

order 

steps 

func 

evals 

backsolves 

LU 

decomps 

accurate 

digits 

time 

Numerical  Results  for 

Problem  1 

GEAR 

io1 

4 

51 

110 

76 

11 

2.9 

2.7 

10  4 

4 

69 

147 

104 

14 

3.6 

3.1 

blended  new  DIFSUB 

103 

4 

33 

98 

134 

10 

3.1 

2.6 

10‘3 

7 

65 

158 

236 

13 

3.8 

5.5 

10 
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10 


-4 


10 

10 
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10 


-3 


10 

10 
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Numerical  Results  for  Problem  2 
GEAR 


4 

96 

211 

130 

20 

2.2 

5.2 

5 

146 

283 

190 

23 

3.4 

8.9 

blended 

new  DIFSUB 

4 

65 

237 

272 

25 

2.7 

7.2 

6 

99 

283 

396 

21 

3.8 

10.0 

Numerical  Results  for  Problem  3 


(4) 


GEAR 

(1001)  (1385)  (1060) 


(54) 


(1.5)  (59.0) 


"too  much  work";  Integration  stopped  at  t  =  8.1 
blended  new  DIFSUB 


5 

7 


139 

210 


382 

520 


570 

822 


16 

18 


2.4 
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Numerical  Results  for  Problem  4 


GEAR 


10-4 

3 

111 

10 

5 

174 

-9 

10 

(5) 

(311) 

-3 

III, 

min 

too  la 

10  J 

6 

166 

257 

358 


164 

233 


23 

31 


1.8 
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blended  new  DIFSUB 
(1376)  (1782)  (121) 


(-1.2) 

;  Integration  stopped  at  t  »  956 
503  716  36  2.5 

Numerical  Results  for  Problem  5 

GEAR 


19.0 

28.0 


6.4 

11.0 


(44.0) 

17.0 
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67 
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129 
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0.1 

5.0 
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240 
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8 
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7.2 
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new  DIFSUB 
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6 

71 

311 

380 

30 

1.6 

10.0 

10  J 

8 

79 

310 

394 

28 

2.9 
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We  had  planned  to  Implement  the  blended  formulas  In  GEAR.  We 
made  changes  to  GEAR  so  as  to  avoid  fixed  dimension  arrays  and  local 
storage.  Improvements  were  made  to  the  calling  sequence  incorporating 
some  ideas  from  the  user  interface  standard  of  Hindmarsh  (1978).  We 
also  added  the  initial  stepsize  selection  algorithm  of  Shampine  and 
Stetter  in  the  source  code  of  RKSW,  which  is  available  as  a  microfiche 
supplement  to  Shampine  and  Wisniewski  (1978).  We  abandoned  this  effort 
when  Hindmarsh' s  new  code  LSODE  became  available.  This  code  already 
has  the  desirable  changes,  and  we  are  putting  the  blended  formulas  into  it. 

A  catalog  of  better  FORTRAN  codes  for  IVPs  has  been  compiled. 

It  has  been  included  in  the  "Working  Papers  of  the  SIGNUM  Meeting  on 
Numerical  ODEs"  and  in  the  June  1979  issue  of  the  SIGNUM  Newsletter. 

A  device  has  been  developed  for  improving  the  convergence  of 
the  corrector  iteration  in  cases  where  coefficient  in  the  decomposed 
matrix  is  out  of  date.  Suppose  we  want  to  solve 


hp'  +  A  -  hf (t,  p  +  BA)  =  0 

for  A  given  a  triangular  factorization  of  I  -  r  ^h8J  where  J  ~  f^  and 
r  is  the  ratio  of  the  current  value  of  h£  to  an  old  value  of  h0.  A 
generalization  of  the  usual  correction  iteration  is  given  by 


A(«fl)  "  A(m)  *  Um[I  ‘  r’lhf3j]_1  residual  (A  (m)) 

where  we  have  introduced  the  relaxation  factor  oj  instead  of  using  a  factor 

m 

of  unity.  For  the  test  equation  f(t,  y)  =  Xy  +  g(t)  with  J  *  X  the 

iteration  error  t.  .  :=  A,  .  -  A  satisfies 
(m)  (m) 

E(mfl) •"  [r  -  y]'ltr(l  -  <V  -  (1  ‘  V)yl£(m) 

where  u  :*  8hX.  If  we  do  not  know  the  number  of  iterations  beforehand, 
we  could  try  to  optimize  the  error  reduction  for  each  iteration  separately 


In  which  case  to  is  the  same  value  to  for  each  iteration.  (If  there  is 
m 

at  least  two  iterations,  we  would  optimize  to^  and  together  and  then 

o>2  and  then  <o^,  etc.  One  might  also  require  100%  error  reduction  for 

X  »  0.)  If  we  restrict  X  so  that  Re  X  £  0  then  Re  y  <_  0  for  the 

methods  of  interest.  The  convergence  factor  is  analytic  for  Re  u  <  0, 

and  so  by  the  maximum  modulus  theorem  the  maximum  value  of  the 

convergence  factor  for  all  Re  h  X  <  0  is  given  by 

sup  )r(l  —  Co)  —  (1  —  wr)iv|/|r  -  ivj  . 

-oo<v<+00 

(It  may  be  more  appropriate  to  compute  the  maximum  over  all  h  in  the 
intersection  of  the  left  half  plane  with  the  absolute  stability  region 
of  the  formula. )  The  square  of  the  worst  case  convergence  factor  is 

max  (r2(l  -  oi)2  +  (1  -  wr)2v2)/(r2  +  v2)  . 


Differentiation  with  respect  to  v  yields  extrema  at  v  =0  and 
2 

v  «  °o.  Therefore  the  worst  case  convergence  factor  is 

max  {  |l  -  (i> | ,  |l  -  rw| } 

for  which 


minimum 


.  Jr  -  1|  „  „  .  _2_  . 

r+1  r  +  1 


This  is  smaller  by  a  factor  of  1/ (r  +  1)  than  the  convergence  factor  for 
u  ■  1.  Experiments  were  performed  with  LS0DE  for  problems  1,  2,  and  3,  but 
the  results  were  inconclusive.  This  device  was  discovered  independently 
by  Chipman  (1979). 

The  Nordsieck  implementation  of  the  blended  formulas  obscures 


the  Identity  of  the  predictor  formulas  that  are  used,  and  interest  has 
been  expressed  in  this  question.  Using  the  ideas  of  Skeel  (1980,  section  2, 
paragraphs  3,  4,  5), we  determined  predictor  formulas  for  the  k-step 
blended  formulas  for  k  ■  1,  2,  3.  The  predictor  for  k  ■  1  is  the  Euler 


formula. 


and  for 


For  k  ■  2  it  is 


-  y  +  y  ,  + 

■'n  n-1 


+  6YhJ{- 


c  =  3  it  is 


23 


-  y  +  y  ,  +  7#  V 

n  n-1  12  n- 


+  ^ ' -  §  yn  + 


,26,  i 

+  15  hVl 


+  64 (yhJ)2{-  ~  yn 


K>|U> 
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4.  Equivalent  Forms  of  Variable  Step  Multlstep  Formulas  (Skeel;  Vu) 

A  manuscript  with  this  title  is  being  prepared  for  publication. 

Widely  available  codes  such  as  A.  C.  Hindmarsh's  LSODE  use  Nordsieck's 

interpolatory  technique  to  vary  the  stepsize.  It  is  often  stated  that 

this  technique  does  not  yield  truly  variable  step  formulas.  However, 

we  have  shown  that  this  is  not  true  in  the  sense  that  there  exists  a 

formula  depending  only  on  the  meshpoints  t  ,  t  . .  t  ,  which  relates 

n  n-1  n— k 

the  computed  values  y  ,  y  y  ,  and  the  derivatives  y' ,  y' 

n  n-1  'n-k  ■'n  n-1 

y^_k»  This  result  may  be  useful  in  improving  the  estimation  of  local 

errors  in  codes  that  use  the  interpolatory  technique.  Other  codes, 

notably  EPISODE, are  based  on  "natural"  variable  step  Adams  and  backward 

differentiation  formulas.  The  implementation  uses  the  scaled  derivatives 

of  the  "modifier  polynomials"  associated  with  these  formulas.  We  were 

interested  in  the  result  of  "blending"  these  two  modifier  polynomials, 

(3) 

for  example,  the  AMF  modifier  polynomial  of  degree  2  at  t  plus 

n 

(2) 

-  h  y  J  times  the  BDF  modifier  polynomial  of  degree  2  at  t  .  After 
n  n  n  n 

a  great  deal  of  algebra  we  determined  the  equivalent  multistep  formula 
satisfied  by  the  computed  solution  and  derivative  values: 

{AMF*3J}  -  h  y  J  {BDF^} 
n  n  n  n  n 


(h  +h  ,)h  . 

n  n-1  n-1 


2 —  <hn  .J  .  ~  h  .Y  J  HI  +  3h  _Y  ,J  1)_1(AMF(2)} 

*■  n-z  n-l  n-l  n-1  n  n  n-Z  n-1  n-1  n-1 


This  is  not  a  blend  of  AMF^  and  BDF^2^  even  for  constant  J  unless  Y 

n 

are  chosen  in  a  manner  incompatible  with  good  stability  behavior. 


n 
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IV.  MULTIGRID  METHODS  (Van  Rosendale;  Skeel) 

A  Ph.D.  thesis  on  this  topic  was  written  by  Van  Rosendale  (1980). 

The  abstract  follows: 

This  thesis  is  concerned  with  the  use  of  multi-level  methods 
to  solve  the  linear  systems  arising  from  finite  element 
discretizations  of  elliptic  equations.  In  all,  three  multi¬ 
level  methods  are  considered.  The  first  of  these  is 
applicable  only  to  quasi-uniform  grids,  but  is  simpler  than 
other  algorithms  considered  in  previous  theoretical  work. 

The  other  two  algorithms  are  applicable  to  both  quasi-uniform 
grids,  and  locally  refined  grids,  those  grids  on  which  the 
size  of  the  largest  and  smallest  elements  may  differ  by  an 
arbitrarily  large  factor.  All  three  algorithms  are 
asymptotically  optimal,  producing  good  solutions  in  0(N) 
operations  on  a  finite  element  grid  with  N  elements.  These 
asymptotically  optimal  complexity  bounds  for  the  last  two 
algorithms  are  the  first  such  bounds  for  multi-level  methods 
on  locally  refined  grids.  One  of  these  algorithms  achieves 
this  0(N)  complexity  bound  under  weaker  than  expected 
conditions  on  the  dimensions  of  the  finite  element  spaces 
used  by  the  algorithm. 

The  multi-level  convergence  results  for  locally  refined  grids 
shown  here  are  based  on  a  new  approximation  result  given  in 
this  thesis.  This  approximation  result  is  of  interest  for 
several  reasons,  the  main  one  being  that  it  is  completely 
local,  making  no  use  of  global  properties  such  as  the 
regularity  of  the  problem.  In  consequence,  it  provides  an 
independent  demonstration  of  the  asymptotically  optimal 
complexity  of  multi-level  algorithms  on  non-convex  domains, 
shown  previously  by  Bank  and  Dupont.  It  also  permits  one 
to  determine  explicit  upper  bounds  on  the  rate  of  convergence 
of  multi-level  methods  on  irregular  finite  element  grids 
using  only  local  properties  of  the  finite  element  space  involved. 

The  table  of  contents  is  as  follows: 

1.  Introduction 

1.1.  Scope. of  Thesis 

1.2.  History  of  Multi-Level  Methods 

1.3.  The  Finite  Element  Approach 

2.  Preliminaries 

2.1.  Elliptic  Equations 

2.2.  Finite  Element  Spaces 

2.3.  Linear  Equations 

3.  Quasi-Uniform  Grids 

3.1.  Introduction 

3.2.  L^  Convergence 

3.3.  Computational  Cost 
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4.  Locally  Refined  Grids 

4.1.  Introduction 

4.2.  Notation 

4.3.  Algorithms 

4.4.  Complexity 

4.5.  Interpolation 

4.6.  Approximation 
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